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coordinate  system  in  which  the  normal  axis  always  coincides  with  the  bisector  of  the  comer.  These  extra¬ 
polation  formulae  are  consistent  with  boundary  conditions  to  first-order  accuracy  in  spatial  increment,  same 
as  the  classical  one-sided  explicit  approximation  scheme  widely  used  for  flat  fiee-surface  case.  Testing 
results  indicate  that  the  present  scheme  works  stably  for  fairly  complicated  geometric  shapes  consisting  of 
ridges  and  valleys  with  steep  and  gentle  slopes  over  a  range  of  Poisson  ratios  of  practical  interest  thus  ena¬ 
bling  us  to  study  more  realistic  problems  for  which  the  topography  plays  a  significant  role  in  shaping  the 
wavcfield  and  analytical  solution  might  not  be  available. 


SECTION  2 

TELES EISMIC  SPECTRAL  AND  TEMPORAL  M0  AND  y.. 

ESTIMATION  FOR  FOUR  FRENCH  EXPLOSIONS  IN  SOUTHERN  SAHARA 

Estimates  for  explosion  moment,  M0,  and  reduced  displacement  potential,  ML.,  are  made  for  four 
French  explosions  at  Taourirt  Tan  Afella  Massif  in  southern  French  Sahara  using  data  from  the  LRSM  net¬ 
work  and  the  arrays  EKA  and_YKA.  Preparatory  to  determining  moments,  /*  estimates  are  made  for  each 
station  and  the  source  region  t*  values  of  0.30  to  0.35  seconds  are  found  for  the  southern  Sahara  test  site. 
This  source  region  attenuation  level  is  consistent  with  the  "hot  spot"  hypothesis  for  the  Ahaggar  plateau  in 
northern  Africa.  Consequently,  the  attenuation  bias  between  Ahaggar  and  the  Nevada  Test  Site  should  be 
small. 

Both  spectral  estimation  and  broadband  temporal  deconvolution  methods  are  used  for  estimation  of 
the  explosion  moments.  The  deconvolution  estimates  of  static  moment  are  found  to  be  consistent  with  the 
spectral  estimation  methods.  Deconvolved  seismograms  for  the  explosions  EMERAUDE,  RUBIS, 
SAPHIR,  and  GRENAT  show  evidence  of  strong  anisotropic  free  surface  interaction  that  may  be  due  to 
scattering  from  the  steep  topography  of  the  Taourirt  Tan  Affela  Massif  test  site. 


SECTION  3 

SCATTERING  FROM  NEAR-SOURCE  TOPOGRAPHY:  TELESEISMIC  OBSERVATION 
AND  NUMERICAL  2-D  EXPLOSIVE  LINE  SOURCE  SIMULATIONS 

2-D  linear  elastic  finite  difference  simulations  of  teleseismic  P  waveforms  from  line  sources  have 
been  used  to  explore  the  variations  that  may  be  induced  in  event  magnitude-yield  determination  by  the 
emplacement  of  explosive  sources  under  mountainous  topographic  features.  The  southern  Sahara  French 
test  site  in  Algeria,  at  Taourirt  Tan  Afella  Massif  on  the  Ahaggar  plateau  has  been  used  as  a  case  study. 
The  topography  of  this  test  site  is  extreme  and  the  event  locations  permit  a  test  of  the  hypothesis  that 
topography  influences  short-period  event  magnitudes  of  contained  nuclear  explosions.  The  maximum  varia¬ 
tion  that  is  expected  is  plus  or  minus  0.15  magnitude  units  from  the  network  mean.  The  magnitude  varia¬ 
tions  are  expected  to  change  rapidly  with  takeoff  angle  and  azimuth. 

Teleseismic  observations  of  the  explosions  at  the  southern  Sahara  test  site  are  compared  to  predic¬ 
tions  made  from  2-D  simulations.  Waveform  data  from  the  arrays  EKA,  and  YKA  ax  well  as  LRSM  data 
have  been  deconvolved  to  broadband  displacement  for  inspection  of  the  apparent  far-field  P-wave  source. 
Qualitative  comparisons  are  favorable  that  the  topography  above  the  explosions  RUBIS  and  SAPHIR 
defocused  teleseismic  pP  at  certain  takeoff  angles  and  azimuths.  Long-period  positive-polarity  pulses  can 
be  seen  at  several  sites  that  may  indicate  Rayleigh-to-P  scattering  from  topography  near  the  source. 

WWSSN  maximum  likelihood  magnitude  data  for  the  "a",  "ab",  and  "max"  P  phases  have  been  used 
to  estimate  that  the  magnitude  variation  due  to  topographic  scattering  is  no  more  than  0.15  rms  magnitude 
units  across  the  WWSSN. 
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SECTION  1 

BOUNDARY  CONDITIONS 
FOR  ARBITRARY  POLYGONAL  TOPOGRAPHY  IN 
A  2-D  ELASTIC  FINITE-DIFFERENCE  SCHEME 
FOR  SEISMOGRAM  GENERATION 

R.-S.  Jih,  K.  L.  McLaughlin,  Z.  A.  Der 
Teledyne  Geotech  Alexandria  Laboratories 
314  Montgomery  Street 
Alexandria,  VA  22314 

INTRODUCTION 

The  2-D  linear  finite-difference  (LFD)  method  has  become  a  popular  means  to  generate 
synthetic  seismograms  because  of  the  ease  with  which  they  can  be  applied  to  model  the  low- 
frequency  response  of  complex  structures  for  which  analytical  solutions  can  not  be  derived.  So 
far  most  applications  were  limited  to  simple  geometric  shapes,  however.  For  instance,  a  lot  of 
earlier  work  only  deal  with  restrictive  models  in  which  the  boundaries  are  either  parallel  to  the 
coordinates  or  only  allowing  45  deg  ramp  (Alterman  and  Karal,  1968;  Alterman  and 
Loewenthal,  1970;  Munasinghe  and  Famell,  1973;  Ilan  et  a!.,  1975,  1979;  Fuyuki  and  Nakano, 
1984;  Hong  and  Bond,  1986;  etc).  Boore  et  al.  (1981)  simulated  vertically  incident  body 
waves  impinging  on  a  45  deg  ramp  by  using  the  same  formulation  as  in  Ilan  et  al.  (1975). 
Ilan  (1977)  utilized  the  LFD  method  to  simulate  the  P-SV  wave  propagation  in  an  elastic 
medium  with  polygonal  free-surface,  but  the  treatment  of  the  transition  points  between  the  seg¬ 
ments  of  various  slopes  in  the  polygons  was  not  addressed.  His  approach  also  requires  the  use 
of  a  non-uniform  grid  and  thus  introduces  some  inevitable  complexity  in  implementing  the 
iteration  procedure.  We  have  found  that  the  treatment  of  the  points  between  the  segments  of 
the  polygon  is  not  a  trivial  problem,  and  that  a  simple  extrapolation  of  the  boundary  conditions 
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from  either  side  of  the  transition  point  as  suggested  in  some  previous  works  may  either  lead  to 
numerical  instability  or  yield  solutions  violating  Huygen’s  principle  and  Snell’s  law. 

In  the  present  work,  we  propose  an  improved  representation  of  boundary  conditions  for 
miscellaneous  types  of  comer  points  on  the  topography  while  maintaining  the  simplicity  of 
fixed  mesh  widths.  For  each  point  on  the  topography  where  slope  changes,  we  use  a  local  grid 
system  in  which  the  normal  axis  always  coincides  with  the  bisector  of  the  comer.  These  extra¬ 
polation  formulae  are  consistent  with  boundary  conditions  to  first-order  accuracy  in  spatial 
increment,  same  as  the  classical  one-sided  explicit  approximation  scheme  for  flat  ffee-surface 
that  has  been  widely  used  for  a  long  time.  Testing  results  indicate  that  the  present  scheme 
works  stably  for  fairly  complicated  geometric  shapes  consisting  of  ridges  and  valleys  with 
steep  and  gentle  slopes  over  a  range  of  Poisson  ratios  of  practical  interest. 


MODEL 


Consider  a  half-space  with  an  arbitrary  polygonal  ffee-surface.  Define  a  Cartesian  coordi¬ 
nate  system  with  X  horizontal  and  Z  vertical,  positive  downward  (Figure  1).  Assume  that  the 
material  is  perfectly  elastic,  isotropic,  and  homogeneous  with  compressional,  shear  velocities 
a,  (3  respectively.  Let  U,  W  be  the  horizontal  and  the  vertical  components  of  displacements. 
Then  the  wave  propagation  is  governed  by  the  following  equation 


d2u(*.z,t)  =  a2  a2u(x,z,t)  +  (az  _  q2  )  a2w(x,z,o  +  ^2  (la) 

dt2  9x2  dxdz  dz2 

a2w(*.z.t)  =  a2  a2w(x,z,t)  +  (a2  _  «2 }  a2u(*.s.t)  +  p2  a^x.z^  (lb) 

at2  9x2  dxdz  az2 

A  grid  is  imposed  on  this  X-Z  plane  with  fixed  mesh  size  Ax  and  Az  .  We  let 


x  =  iAx,  z  =  kAz  and  t  =  LAt  where  At  is  the  increment  step  in  time,  i,k,l  are  positive  integers, 
and  denote  by  Ul  kj  and  Wjkj  the  approximate  components  of  displacement  at  grid  point 
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(iAx.kAz)  at  time  lAt.  By  replacing  the  derivatives  in  equation  (1)  by  central  finite  differences, 
we  obtain  the  following  explicit  formulae: 

2 

Uijtj+i  =  -  Ujjcj.  i  +  2Uiikii  +  a2(-^-)  (  Uj+ijcj  -  2Uitkj  +  Ui.iw  ) 

2 

+  P2(-|^)  (  Ui>k+1J  -  2Ui>k>1  +  UiJc_u  ) 

+  a2  -  P2) .  (  Wi+U+U  -  Wi+U_1(!  -  WiW  Ww>wj)  (lc) 

2 

Wiik4+,  =  -  WU4_,  +  2Wltk(]+p2(^-)  (  Wj+UJ  -  2WiJcJ  +  Wi_lw  ) 

+  a2(-^H  W,k+U  -  2WW  +  Wi(k_J(1 ) 

+  «2  -  P2)  •  (  ui+t.k+l.l  -  Uwjc-14  -  Ui_u+u+  Uw.fc.jj)  (Id) 

The  reader  is  referred  to  Boore  (1972),  Alterman  et  al.  (1972),  and  Kelly  et  al.  (1976)  for 
a  more  detailed  discussion  of  the  approximation  of  partial  derivatives  by  finite  differences. 
Recently  a  number  of  techniques  have  been  pursued  in  an  effort  to  improve  the  computational 
performance  of  LFD  solutions  to  wave  equations.  These  include  variable  grids,  higher  order 
schemes  or  implicit  rather  than  explicit  methods.  None  of  these  can  provide  significant 
improvement  over  the  traditional  explicit  second  order  centered  LFD  schemes  described  by 
Alterman,  Boore  and  Kelly  when  complicated  topographic  profiles  are  used  as  ffee-surface.  To 
implement  arbitrary  topography  in  the  LFD  scheme,  it  suffices  to  approximate  the  topography 
by  polygons.  Six  algorithms  are  presented  here  to  implement  the  free-surface  boundary  condi¬ 
tion  for  the  six  separate  cases  shown  in  Figure  1  labeled  (A)  through  (F).  Only  the  formulae 
for  increasing  elevation  with  increasing  x  are  given  for  brevity.  Since  the  calculation  of  comer 
points  always  requires  displacements  of  points  nearby,  it  is  necessary  to  follow  a  specific  order 
in  each  iteration  step  once  the  displacements  at  flat  free-surface  are  extrapolated.  This  specific 
order  is  suggested  in  Figure  1  and  self-explanatory  from  the  description  of  various  algorithms 
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in  the  following  sections.  Also  the  topography  must  be  sampled  so  that  each  segment  of  the 
polygon  consists  of  at  least  3  points  (i.e.  2  sub-segments). 


ALGORITHMS 


(A)  CONSTANT  STEEP  SLOPE 

Suppose  that  a  grid  point  <i,k>  lies  on  an  inclined  free-surface  with  slope 
mAzJAx,  m  >  2.  The  boundary  condition  states  that  both  the  tangential  and  normal  stress 
components  vanish  at  grid  point  <i,k>,  i.e.  if  we  rotate  the  X-Z  coordinate  system  by  angle 
0  =  tan~*(m  Az/Ax)  counterclockwise  to  the  X-Z'  system,  then 


(2a) 

(2b) 


U'ije  =  U'0  -  cos0sine(WVU-,-W'i(k+m-i)  (3a) 

W'i>k  =  W'0  -  cos0sin0(  l-2-^)(U/i+1,k_,-U'i>k+m_1)  (3b) 

where  U0  =  cos20  Ui  k+m  +  sin20Ui+i>k,  W0  =  cos20  WiJc+m  +  sin20Wi+1  k  are  the  linearly 
interpolated  displacements  at  the  point  labeled  "0"  (Figure  2).  Rotate  the  displacement  at 
points  "0",  <i+l,k-l>,  and  <i,k+m-l>  to  the  inclined  system  to  obtain  U'ik,  W'jk  by  using 
formulae  (3).  The  displacements  at  <i,k>  are  then  rotated  back  to  the  original  system  by  angle 
-0  . 
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(B)  CONSTANT  GENTLE  SLOPE 

For  grid  point  <i,k>  on  an  inclined  free-surface  with  slope  AzJAx  (Figure  3),  let 

L0  =  EoUi+1>k  +  (  1  -  Eo  )Uitk+,  (4a) 

W0  =  eoWi+1)k  +  (  1  -  )Wi>k+1  (4b) 

where 


Compute 


Eo  =  sin20  , 

0  =  tan-1(Az/Ax), 


U'ilk  =  U'(j-sin0cos0(W'i+i.k  -  W'ik+1)  (5a) 

W'y.  =  W'o-sin0cos0(l  -  2-^)(U'i+1Jt  -  U',  k+1)  (5b) 

or 

Then  rotate  U'j  k,  W'ik  back  to  the  original  system.  Note  that  this  algorithm  is  simpler 
than  Ilan’s  (1977)  DuFort-Frankel  type  approximation  (see  Appendix)  and  the  weighting  factor 
sin0cos0  in  formulae  (5)  is  applicable  whether  0  S  45°  or  not.  Also  note  that  the  words 
’gentle’  and  ’steep’  have  only  relative  meaning  here  since  we  can  always  adjust  Ax  and  Az 
before  implementation. 


(C)  CONCAVE  HORIZONTAL-TO-GENTLE  SLOPE  TRANSITION 

For  grid  points  at  the  bottom  of  a  canyon  with  slope  mAz/Ax  ,  m  £  1  (Figure  4),  assum¬ 
ing  that  Ax  ^  Az  ,  let  0  =  tan_1(mAz/Ax),  Eq  =  tan(0/2)Az/Ax,  Ei  =  tan(0/2)cot0,  and 


Uo  -  (l-£o)Ujjk+i  +  £<)Uj+ijk+i  (6a) 

W0  =  (l-e0)Wlk+1 +e0W1+Ik+I  (6b) 

Ui  =  (l~Ei)U,_ltk+m+1  +  £|Uj_ik+1  (6c) 

wi  =  (l-eOWj.i  k+m+i  +  £)  Wj_|  k+|  (6d) 
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U2  =  (l-E^Ui+ijc-nH.!  +  EiUi+1Jt+1 

(6e) 

w2  =  (l-£i)Wi+1>k_IIH.1  +  E]Wi+1Jc+1 

(6f) 

Then, 

u'i'k=u'°-|:(W'rW',) 

(7a) 

W'u  =  W'o  -  ^-(i  -  2-^kuvu'O 

2Ax  a 1 

(7b) 

Note  that  we  rotate  U,  W  by  angle  0/2  to  get  U',  W'  (instead  of 

0  )  so  that  Z'  axis  is 

consistent  with  the  direction  of  the  local  bisector  at  grid  point  <i,k>. 

For  the  case  Ax  <  Az  ,  substitute  Eq  =  cot(0/2)Ax/Az,  and  £j  =  tan0  tan(0/2)  in  formu¬ 
lae  (6),  and  replace  Az/2Ax  in  (7)  by  cot(0/2)/2  . 


(D)  CONCAVE  GENTLE-TO-STEEP  SLOPE  TRANSITION 

For  grid  points  with  left  slope  nAz/Ax  and  right  slope  mAz/Ax,  m  >  n  (Figure  5), 
assuming  that  Ax  £  nAz  ,  let  0  =  tan-I(nAz/Ax),  <t>  =  tan_1(mAz/Ax),  T|  =  (0  +  4>  )/2, 
Eo  =  Ei  =  tanOcott),  e2  =  tanT\cot0,  and 


Then, 


Uo  =  (l-Eo)Ui+iilt  +  EoUj+ijc+n 

(8a) 

W0  =  (l-Eo)Wi+l,k  +  ^Wi+ljc+n 

(8b) 

Ui  =  (l-Ei)Uijk+n  +  EjUj-iij+n 

(8c) 

Wj  =  (l-Ei)Wi_llc+n  +  £iWj_ljk+n 

(8d) 

U2  =  (l-Ei)Ui+lk  +  EiUi+1>k_m 

(8e) 

W2  =  (l_el)^i+l,k  +  el^i+l,k-m 

(80 

U'u  =  U  o  - 


1 


tanr|+tan0 


(W'2-W',) 


(9a) 
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W'u  =  W'0 - - - -(1  -  2^-r )(U'2-U'i) 

0  tanrj+tan0  a2 


(9b) 


Note  that  we  rotate  U,  W  by  angle  Tj  to  get  U',  W'  so  that  the  Z'  axis  is  consistent  with 
the  direction  of  the  local  bisector  at  grid  point  <i,k>. 


For  the  case  that  Ax  >  nAz,  substitute  £q  =  tanGtanr),  e,  =  tan0cotr|,  and  e2  =  tant|cot<li, 
in  formulae  (8),  and  replace  tarn)  +  tan0  in  (9)  by  cott|  +  cot0. 


(E)  CONVEX  CHANGE  IN  SLOPE 

For  grid  points  with  change  in  slope  from  mAz/Ax  to  nAz/Ax,  m  >  n  ,  let 
0  =  tan~'(nAz/Ax),  <j>  =  tan-1(mAzMx),  r|  =  (0  +  $  )/2,  Eq  =  tantjAz/Ax.  If  Eq  is  less  than  1, 
then  let 


otherwise  let  e3  =  I/Eq  ,  and 


Uo  =  (l-£o)Ui.k+l  +  ^O^i+l.k+1  (10a) 

W0  =  (l-€o)WjJc+1  +  EoWl+lk+1,  (10b) 

Uq  =  (l-C3)Ui+l,k+l  +  e3Ui+l,k  (10c) 

W0  =  (i-e3)Wi+ik+j  +  EsW^.,  k  (lOd) 


Now  depending  on  the  value  of  m  ,  we  have  2  cases  : 

(E.l)  m=l  (Figure  6A),  note  that  tan<j>  <  tanr)  <  tan0  ,  hence 
Ej  =  (tanr)cot<5)  —  1  )/(n — 1 )  is  always  between  0  and  1.  Let 


U,  =  (l-ei)Ujk+1  +  EiUiik+n  (11a) 

W,  =  (l-Ei)Wj  k+1  +  E]  Wj  k+n  (lib) 


Then, 


U'i,k  =  U'0  -  e2(WVilk-W'1) 


(12a) 
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W'u  =  W'0  -  62(1  -  2-^XUV^-U',)  (12b) 

or 

where  c2  =  tan<>  if  Co  S  1  ,  or  cotT)  if  Eo  >  1  . 

(E.2)  For  the  case  m  >  1  (Figure  6B),  let  £j  *  (2Azcotri/Ax) ,  then  Ej  is  always 

between  0  and  1. 

U,  =  (l-£i)Ui,k-i  +  EiUj+i.k-i  (He) 

Wi  =  (l-Ei)Wi.k-l  +  £lWi+Uc-l  (Hd) 

The  displacements  at  point  <i,k-l>  are  computed  by  extrapolation  (see  algorithm  (G))  since  it’s 

outside  the  grid.  Now, 

U'itk  =  U'0  -  e2(W',-W'i>k+1)  (12c) 

W'1>k  =  W'0  -  £2(1  -  2-f£)(U',-U'i.k+1)  (1 2d) 

or 

where  e2  =  0.5tanr|  if  Eq  <,  1  ,  or  Ax/(2Az)  if  Eq  >  1  . 

Note  that  again  we  rotate  U,  W  by  angle  t)  to  get  U\  W'  so  that  the  Z'  axis  is  con¬ 
sistent  with  the  direction  of  the  local  bisector  at  grid  point  <i,k>. 


(F)  CONVEX  GENTLE  SLOPE-TO-HORIZONTAL  TRANSITION 

Without  loss  of  generality,  we  may  assume  that  the  topography  changes  slope  from 
Az/Ax  to  0  at  the  crests  (Figure  7).  Let  $  =  tan_l(Az/Ax),  Co  =  tan($/2)tan<t>,  and 
E[  =  tan(<{>/2)cot<j>.  Note  that  £<3  <,  1  provided  Ay  <,  V3  Ax  ,  i.e.  $  £  60°.  Let 


Co  =  (l-EoJU^k+i  +  EoU^ijt+1 

(13a) 

w0  =  (l-Eo)Wj>k+,  +  EoWi+I>k+, 

(13b) 

U,  =  (l-€i)D|-l,k+l  +  elO,_i,k+2 

(13c) 

W,  =  (1-E,)WMJC+1  +  elWi-l,k+2 

(13d) 
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U2  =  (l-eiHW+1  +  ejUj+ijc  (13e) 

W2  =  (l-€1)Wi+u+1+e,Wi+u  (13f) 

Then, 

U'i(k  =  U'o  -  -^-(W'2-W',)  (14a) 

W'i<k  =  W'0  -  -i2Di(i  -  2-^-KU'2-U',)  (14b) 

2  or 


Here  we  rotate  U,  W  by  angle  <J>/2  to  get  U',  W'  so  that  the  Z'  axis  is  again  con¬ 
sistent  with  the  direction  of  the  local  bisector  at  grid  point  <i,k>. 

(G)  FICTITIOUS  BOUNDARY  POINTS  FOR  STEEP  SLOPE 

Suppose  that  a  grid  point  <i,k>  lies  on  the  ffee-surface  with  slope  mAz/Ax,  1  .  Due 
to  the  explicit  scheme  (equation  (lc)  and  (Id))  we  adopt  for  the  interior  points,  the  calculation 
of  the  displacements  at  inner  points  <i+l,k>,  <d+l,k-l>,...,<i+l,k-m+l>  require  displacement 
values  at  <i,k-l>,  <i,k-2>,...<i,k-m>  which  are  outside  the  topography.  This  difficulty  can  be 
easily  removed  by  evaluating  m  fictitious  layers  above  the  topography  as  follows  (Figure  8): 

For  1  £  n  £  m,  let  Eq  =  (n/m)sin29  ,  where  9  =  tan-1  mAz/Ax  ,  use 

U0  =  eoU1+lik_m  +  (l-eo)Ui>k  and  W0  =  EoWi+lk_m  +  (l-eo)Wik  as  the  approximate  displace¬ 
ment  of  the  orthogonal  projection  of  point  <i,k-m+n>  on  the  ffee-surface. 

Now  rotate  displacements  at  points  <i+l,k-m>,  <i,k>  and  "0"  by  angle  9  as  before  to 
extrapolate  the  approximate  motion  at  point  <i,k-m+n>  : 

U'i.k-™  =  UV—  sin9cos9(W'I+1(k_m  -  W'u) 
m 

WWn  =  W'0  -  (1-2-^-)— sin9cos9(U'i+l  k_m  -  U'j>k) 
a  tn 
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Then  rotate  U',  W'  back  to  the  original  system. 


EXAMPLES 

EXAMPLE  1 

Figure  9  shows  the  propagation  of  a  normally  incident  P  wave  through  a  model  with  a 
45°  ramp  on  the  top  of  grid.  Shading  is  proportional  to  displacement  amplitude.  The  initial 
wave  was  a  tapered  displacement  pulse  of  the  form  te-011  and  von  Neumann  (i.e.  0-slope  or 
symmetric)  boundary  condition  was  used  on  both  sides.  The  appropriate  P-S  conversions  and 
the  reflections,  diffractions  satisfying  Snell’s  law  and  Huygen’s  principle  are  clearly  visible  in 
these  successive  snapshots  taken  every  sec. 

EXAMPLE  2 

Figure  10  shows  the  snapshots  of  displacement  fields  generated  by  an  upgoing  P  wave  of 
velocity  5.2  km/sec  in  a  grid  of  size  6.5km  by  3.9km  with  steep  topographic  configuration. 
The  successive  frames  separated  by  0.125  sec  show  the  initialization  of  the  wave  and  locations 
of  sensor  1  to  7  (left  array)  ,  8  to  17  (right  array)  which  were  separated  50  meters  apart  verti¬ 
cally  (A),  P-reflection  followed  by  S  wave  starting  at  right  (B,C),  completely  developed 
reflections  from  all  parts  of  the  topography  (E)  and  complex  wavefields  containing  reflections, 
diffractions  and  possibly  excited  surface  waves  (E,F).  The  initial  P  wave  has  an  incident  angle 
of  20°  and  the  topography  is  a  (due  north  344°)  cross-section  of  Taourirt  Tan  Afella  Massif  in 
southern  Algeria  (Faure,  1972).  This  granite  mountain  rises  from  the  (1100  meter  elevation) 
Ahaggar  plateau  to  a  peak  of  2100  meters  in  less  than  1  km  and  has  slopes  that  locally  exceed 
53°.  Several  nuclear  tests  were  conducted  here  from  1962  to  1966.  To  cover  a  wide  range  of 
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slope  in  LFD  calculation,  a  grid  with  horizontal  spacing  (Ax)  0.05  km  and  vertical  spacing 
(  Az )  0.01  km  was  used.  Figure  11  shows  the  dilatational  strains  recorded  at  17  sensors 
buried  beneath  the  peaks  of  the  ridges.  It  can  be  observed  that  the  free-surface  reflection  is 
severely  altered  due  to  scattering  from  the  free-surface.  The  symmetric  boundary  conditions  on 
the  sides  also  causes  some  undesirable  edge  reflections  from  left  side.  For  detailed  discussion 
of  this  example  and  its  applications  in  modeling  teleseismic  P  waveforms  observed  at  LRSM 
stations,  see  section  3  of  this  report. 

EXAMPLE  3 

Figure  12  and  13  show  a  Rayleigh  wave  incident  on  an  upgoing  and  a  downgoing  ramp 
of  height  h  -  6.4  km  =  64  Az  respectively.  The  initial  waveform  is  chosen  so  that  the  vertical 
displacement  W  on  the  flat  free-surface  is  a  Ricker  wavelet.  This  wavelet  has  been  adopted 
frequently  in  finite-difference  simulations  because  it  is  well  localized  in  both  spatial  and 
wavenumber  domains  (Boore,  1972;  Munasinghe  and  Farnell,  1973;  Fuyuki  and  Matsumoto, 
1980;  Fuyuki  and  Nakano,  1984;  Levander,  1985).  To  avoid  the  grid  dispersion,  we  chose  3.2 
km  =  32  Ax  as  the  dominant  wavelength  of  the  incident  Rayleigh  wave  packet  in  a  homogene¬ 
ous  medium  of  P  wave  velocity  5  km/sec  and  Poisson  ratio  0.35.  Absorbing  boundary  condi¬ 
tions  given  by  Clayton  and  Engquist  (1977)  and  Emerman  and  Stephen  (1983)  were  used  on 
the  sides  and  bottom  to  suppress  the  artificial  reflections  from  the  sides  of  the  grid.  Figures 
(A)  through  (F)  correspond  to  displacement  wavefields  at  distinct  times  with  a  temporal  spac¬ 
ing  of  2  sec.  In  both  cases  the  comer  points  always  act  as  point  scatterers  radiating  converted 
body  waves  away  in  frames  (C)  through  (F).  Note  that  the  diffracted  pattern  in  the  downgoing 
ramp  case  (Figure  13)  is  much  more  complicated  than  in  the  upgoing  case  (Figure  12)  and  that 
the  quasi-transparent  boundary  conditions  allow  the  wave  to  disappear  into  the  sides  and  bot¬ 
tom  of  the  grid.  Even  a  moderate  topography  can  substantially  attenuate  a  short-period 
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fundamental  Rayleigh  wave  due  to  body  wave  conversion  (McLaughlin  and  Jih,  1986b).  Such 
attenuation  of  short-period  fundamental  Rayleigh  wave  incident  upon  topography  is  usually 
stronger  than  that  due  to  shallow  heterogeneous  sediment  over  homogeneous  half  space 
(McLaughlin  and  Jih,  1987). 

CONCLUSION 

In  this  paper,  we  have  presented  a  straightforward  method  to  simulate  the  P-SV  wave 
propagation  in  2-D  LFD  models  with  free-surface  topography.  This  method  specifies  the  treat¬ 
ment  of  a  polygonal  boundary  including  the  transition  points  between  the  straight  line  seg¬ 
ments.  Boundary  conditions  are  implemented  in  a  way  that  at  transition  points  the  normal  axis 
always  coincides  with  the  bisector  of  the  comer.  This  approach  reflects  the  natural  idea  to 
account  for  the  slope  change  along  the  topography.  It  is  empirically  stable  for  those  Poisson 
ratios  of  interest  even  in  models  with  quite  complex  topography  and  is  suitable  for  the  calcula¬ 
tion  of  teleseismic  P  waves  from  explosive  line  source  under  realistic  surfaces.  The  numerical 
scheme  in  the  present  work  has  truncation  error  of  order  one  in  spatial  increment,  consistent 
with  the  classical  one-sided  explicit  extrapolation  formulae  widely  used  in  the  flat  free-surface 
case.  Schemes  with  higher  order  of  accuracy  might  be  harder  to  derive  because  of  the  asym¬ 
metrical  differences  and  interpolations  involved. 
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FIGURE  CAPTIONS 

FIGURE  1.  Cartesian  coordinate  system  in  2-D  finite-difference  scheme  and  miscellaneous 
comer  points  on  the  topography.  Wave  propagation  is  constrained  to  lie  in  this  X-Z  plane. 
Labels  of  the  points  correspond  to  the  specific  algorithm  to  use.  It  suffices  to  use  combination 
of  such  transition  points  and  points  on  inclined  surface  to  approximate  any  polygonal  topogra¬ 
phy. 

FIGURE  2.  Ramp  with  slope  2Az/Ax.  Same  algorithm  applies  to  any  ramp  with  slope 
mAz/Ax,  m£2. 

FIGURE  3.  Ramp  with  slope  Az/Ax. 

FIGURE  4.  The  bottom  of  a  canyon  of  slope  2Az/Ax.  Same  algorithm  applies  to  any 
horizontal-to-inclined  transition  of  slope  mAz/Ax  where  m  £  1. 

FIGURE  5.  Comer  where  slope  changes  from  nAz/Ax  to  mAz/Az  ,  m  >  n  . 

FIGURE  6A.  Comer  where  slope  changes  from  nAz/Ax  to  Az/Az, 

FIGURE  6B.  Comer  where  slope  changes  from  nAz/Ax  to  mAz/Az  ,  1<  m  <  n  . 

FIGURE  7.  Inclined-to-horizontal  transition. 

FIGURE  8.  Example  of  steep  slope  mAz/Ax,  m  =  4.  The  inner  points  <i+l,k>  to  <i+l,k-4> 
require  fictitious  points  at  <i,k-l>,...<i,k-4>  for  solution  of  equation  (1). 

FIGURE  9.  Snapshots  of  displacement  fields  generated  by  a  normally  incident  P  wave  pro¬ 
pagating  across  a  homogeneous  medium  of  size  39km  by  26km  (390  x  260  mesh  points)  with 
a  45°  ramp  on  the  top.  Frames  are  separated  by  1  sec  and  the  darkness  in  each  frame  is  pro¬ 
portional  to  the  displacement  Note  the  reflection  at  the  inclined  surface  breaks  up  into  P-  and 
S-waves,  and  the  lower  comer  radiates  body  waves  away  circularly  (C). 

FIGURE  10.  Snapshots  of  displacement  fields  with  temporal  spacing  0.125  sec  generated  by  a 
broadband  P  waves  incident  on  a  6.5km  by  3.9km  (130  x  390  mesh  points)  cross-section  of 
Taourirt  Tan  Afella  Massif  with  P  wave  velocity  5.2  km/sec.  The  incident  P  wave  has  an 
incident  angle  of  20°. 

FIGURE  11.  Dilatational  strain  recorded  at  17  sensors  buried  beneath  the  peaks  of  the  topo¬ 
graphic  profile.  Sensors  are  50  meters  apart  vertically  with  locations  shown  by  solid  circles  in 
Figure  10A.  Number  1  through  7  (from  top  to  bottom)  are  on  the  left  array. 

FIGURE  12.  Rayleigh  wave  incident  on  a  6.4  km-high  upgoing  ramp.  Snapshots  show  the 
displacement  wavefields  with  temporal  spacing  2  sec.  The  body  waves  excited  by  Rayleigh 
waves  on  encountering  the  comers  are  clearly  seen  propagating  away  circularly  in  frames  (C) 
through  (F).  Note  that  the  diffracted  P  and  S  waves  propagate  mainly  to  the  right  while  for 
downgoing  ramp  (Figure  13)  the  diffracted  body  waves  propagate  more  isotropically. 

FIGURE  13.  Same  as  Figure  12  except  for  opposite  direction  of  the  ramp.  The  diffracted 
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wavefield  is  more  complicated  as  compared  to  the  case  of  upgoing  cliff  with  the  same  height 
(Figure  12).  Note  the  constructive  interference  of  the  forward-scattered  S  wave  by  the  two 
corners  and  destructive  interference  for  the  back-scattered  S  wave. 

FIGURE  14.  Ilan’s  approach  for  inclined  surface.  A  nonuniform  grid  with  varying  Ax  is 
required  in  order  to  apply  this  algorithm  to  realistic  polygonal  free-surface. 
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APPENDIX 


ILAN’S  METHOD  FOR  CONSTANT  SLOPE  RAMP 


For  grid  point  <i,k>  on  an  inclined  free-surface  with  slope  Az/Ax  (Figure  17),  let 


U'i.k  =U'o-sin20(W'i+1>k  -  W'u+t) 

W'i>k  =W'0-sin28(l  -  ~  U'yc*,) 


where 


Uo  -  +  (  1  -  Eo  )Uj(k+2 

Wq  =  £oWi+lik+1  +  (  1  -  Eq  )Wi(k+2 
eo  =  2sin28  if  8  £  45°, 

8  =  tan-1  (Az/Ax), 

if  8  £  45°,  then  take  Eo  =  2cosz8,  and  replace  equations  (16a)  and  (16b)  by 


Uo  =  eoUi+lik+1  +  (  1  -  Eq  )Ui+2,k  (16c) 

W0  =  EoWj^jjj+j  +  (  1  -  Eo  )Wi+2jc  (16d) 

Then  rotate  U'ik,  W'j  k  back  to  the  original  system.  We  have  observed  that  without 

appropriate  manipulation  of  the  transition  points  on  the  topography  as  proposed  in  the  present 
work,  this  algorithm  alone  would  suffer  from  instabilities. 
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SECTION  2 

TELESEISMIC  SPECTRAL  AND  TEMPORAL 
M0 AND 

ESTIMATES  FOR  FOUR  FRENCH 
EXPLOSIONS  IN  SOUTHERN  SAHARA 

K.  L.  McLaughlin,  A.  C.  Lees,  and  Z.  A.  Der 
Teledyne  Geotech  Alexandria  Laboratories 
314  Montgomery  St. 

Alexandria,  VA  22314 

INTRODUCTION 

Several  contained  nuclear  explosive  tests  were  conducted  in  southern  Sahara  from  1962 
to  1966.  These  tests  are  reported  in  Duclaux  and  Michaud  (1970)  and  Faure  (1972). 
Because  U.S.  testing  experience  in  granite  is  limited,  these  tests  present  an  opportunity  to 
extend  the  U.S.  knowledge  of  tests  in  granite.  In  this  paper  we  report  on  estimates  of  the 
explosion  source  from  short  period  Long  Range  Seismic  Measurements  (LRSM)  stations  and 
the  AWRE  U.K.  arrays  EKA  and  YKA.  We  estimate  the  explosion  source  strength  parameter¬ 
ized  as  the  low  frequency  limit  of  the  seismic  moment,  M0,  or  the  low  frequency  limit  of  the 
reduced  displacement  potential  (RDP),  VF„. 

Schock  et  al.  (1972)  and  Heuze  (1983)  report  on  geophysical  parameters  of  the  Taourirt 
Tan  Afella  Massif  located  in  southern  Algeria.  We  have  adopted  the  values  of  P-wave  velo¬ 
city,  a  -  5.8  km/s,  S-wave  velocity  f3  -  3.3  km/s,  and  density  p  -  2.6  gm/cc  as  representative 
of  the  material  in  the  vicinity  of  the  tests.  The  Taourirt  Tan  Afella  Massif  is  located  on  the 
Ahaggar  plateau  identified  by  Crough  (1981)  as  a  1000  km  wide  dome  of  continental  crust 
uplifted  by  a  mantle  hot  spot.  The  dome  has  an  average  elevation  of  1100  meters  and  is 
capped  by  late  Tertiary  alkali  basalt  volcanos.  The  Taourirt  Tan  Afella  Massif  is  a  16  km2 
granite  mountain  that  rises  to  a  2100  meter  peak  in  a  lateral  distance  of  less  than  1  km.  Given 
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this  tectonic  setting,  the  source  region  t*  estimate  of  0.35  sec  given  by  Der  et  al.  (1985a  and 
b)  is  in  accord  with  an  upper  mantle  of  high  intrinsic  attenuation. 

Figure  1  shows  the  location  of  the  test  site,  several  LRSM  stations  and  the  arrays  EKA 
and  YKA  analyzed  in  this  paper.  Table  1  gives  the  origin  times  and  locations  given  by 
Duclaux  and  Michaud  (1970)  for  the  four  explosions  analyzed  in  this  paper.  Figure  2  shows 
these  locations  superimposed  on  a  topographic  map  from  Faure  (1972). 


TABLE  1. 


DATE 

ORIGIN  TIME 

LAT  (N) 

LONG  (E) 

EVENT 

1963/03/18 

10:02:00.351 

24.0414 

5.0522 

EMERAUDE 

1963/10/20 

13:00:00.011 

24.0355 

5.0386 

RUBIS 

1965/02/27 

11:30:00.039 

24.0587 

5.0312 

SAPHIR 

1966/02/16 

11:00:00.035 

24.0442 

5.0412 

GRENAT 

t*  ESTIMATES  FOR  THE  SOUTHERN  SAHARA  TEST  SITE 

Before  explosion  size  estimates  may  be  made,  the  intrinsic  attenuation  for  the  paths  must 
be  determined.  We  follow  a  procedure  much  like  that  in  Der  et  al.  (1985)  to  partition  the  t* 
estimate  into  a  source  t*  and  a  receiver  t*.  In  order  to  estimate  intrinsic  attenuation  for  the 
available  paths,  a  procedure  was  followed  similar  to  one  described  in  Der  and  Lees  (1985). 
The  seismic  spectra  were  corrected  for  instrument  response  and  approximate  von  Seggern  and 
Blandford  (1972)  RDP.  The  far-field  P-wave  spectra  is  assumed  to  be  proportional  to 

(1+AZK2)1'2 

(l+K2)3'2 

where  A=5.08,  K=tu/k,  and  k=16.8(5/Y)13.  Marshall  et  al.  (1979)  give  the  yields  for  RUBIS 
and  SAPHIR  as  Y=52  and  Y-120  KT  respectively.  The  corrected  amplitude  spectra  in  the 
bandwidth  between  0.5  and  3.0  Hz  were  fit  to  a  form  of  fIoexp(-rt  t*  0-  Signal-to-noise  esti¬ 
mates  for  each  frequency  were  estimated  and  only  those  values  with  an  estimated 
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(signal+noise)-to-noise  power  ratio  greater  than  2  were  used.  The  signal  spectra  were  generally 
absent  of  any  useful  information  for  frequencies  above  3  Hz.  Whole  path  t*  estimates  are  listed 
for  each  station  in  Table  2  along  with  estimated  standard  errors.  Only  RUBIS  and  SAPHIR 
were  large  enough  to  permit  useful  estimates  of  t5.  Figure  3  shows  examples  of  spectra 
corrected  for  instrument  response  and  von  Seggem  and  Blandford  (1972)  RDP. 

We  assume  that  the  whole  path  t*  is  the  sum  of  the  station  t*  plus  the  test  site  t*,  and  we 
wish  to  estimate  the  test  site  t*.  Der  et  at.  (1985a)  list  t*  contributions  to  be  expected  for 
RK-ON,  EKA,  YKA,  HN-ME,  NP-NT,  CPO,  NAO,  and  GRFO.  The  sites  RK-ON,  SB-GR, 
and  OO-NW  are  now  occupied  by  RSON,  GRFO,  and  NAO  respectively.  Given  these  station 
estimates  for  t*p,  we  estimate  t*P  for  Ahaggar  as  0.31  (0.04)  sec.  Der  et  ai  (1978)  list  esti¬ 
mates  for  short  period  S-wave  t*^  for  the  LRSM  stations  BL-WV,  RK-ON,  RY-ND,  VO-IO, 
DH-NY,  HL-1D,  KN-UT,  MN-NV,  FR-MA,  LC-NM,  RT-NM,  and  GI-MA.  Assuming  t*i  = 
4  t*P,  the  test  site  estimate  for  Ahaggar  t*P  is  0.26  (0.05)  sec.  Because  this  value  depends  on 
more  assumptions  relating  t*s  10  t*P,  it  is  probably  less  reliable.  All  these  data  are  consistent 
with  the  Der  et  al.  (1985b)  estimate  of  0.35  sec  for  the  southern  Algeria  test  site. 
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TABLE  2. 


t*  ESTIMATES  FROM  RUBIS  AND  SAPHIR 

STATION 

F 

o 

sec 

sec 

BL-WV 

0.37 

0.04 

BR-PA 

0.53 

0.06 

CPO 

0.53 

0.06 

DH-NY 

0.46 

0.12 

EB-MT 

0.62 

0.04 

EKA 

0.52 

0.10 

EU-AL 

0.44 

0.06 

FR-MA 

0.50 

0.06 

SB-GR 

0.55 

0.06 

GI-MA 

0.48 

0.07 

HL-ID 

0.51 

0.04 

HN-ME 

0.34 

0.04 

KN-UT 

0.69 

0.04 

LC-NM 

0.59 

0.03 

NP-NT 

0.49 

0.07 

OONW 

0.46 

0.06 

RK-ON 

0.37 

0.05 

RT-NM 

0.70 

0.04 

RY-ND 

0.57 

0.06 

VO-IO 

0.49 

0.07 

WO-AZ 

0.80 

0.05 

YKA 

0.45 

0.10 

SPECTRAL  ESTIMATES  FOR  M0,  and 

In  addition  to  the  slope  of  each  spectra,  the  zero  frequency  limit,  12^,,  was  estimated  for 
each  spectra.  These  values  are  tabulated  in  Table  3  along  with  estimates  for  M0  and 
derived  from  geometrical  spreading  given  by  the  J.B.  travel  time  tables  (Aki  and  Richards, 
1980).  In  the  calculations,  surface  values  at  each  station  of  a,  P,  and  p  were  assumed  to  be  a 
=  5.8  km/sec,  (3  -  3.3  km/sec,  and  p  -  2.6  gm/cc,  except  for  the  stations  RK-ON,  EU-AL, 
BL-WV,  KN-UT,  HN-ME,  DH-NY,  HL-ID,  NP- NT,  and  BR-PA  where  (a,  p,  p)  values  were 
obtained  from  Der  et  al  (1977).  These  values  are  also  listed  in  Table  3. 
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TABLE  3. 


RUBIS  AND  SAPHIR  LRSM  SPECTRAL  MOMENT  ESTIMATES 

STATION 

A 

t* 

Ac 

P 

a 

P 

M0 

4^ 

deg 

s 

nm-s 

gm/cc 

km/s 

km/s 

dyne-cm 

m3 

RUBIS 

BL-WV 

72.6 

0.37 

50. 

2.7 

6.0 

3.46 

3.77E+22 

3433. 

DH-NY 

66.7 

0.46 

28. 

2.7 

6.0 

3.46 

2.07E+22 

1889. 

EB-MT 

78.1 

0.62 

102 

2.6 

5.8 

3.3 

8.04E+22 

7321. 

EU-AL 

79.4 

0.44 

189. 

1.9 

2.7 

1.4 

8.93E+22 

8121. 

FR-MA 

86.3 

0.50 

118. 

2.6 

5.8 

3.3 

l.OOE+23 

9136. 

GI-MA 

84.4 

0.48 

97. 

2.6 

5.8 

3.3 

8. 14E+22 

7409. 

HL-ID 

92.2 

0.51 

45. 

2.5 

5.2 

2.96 

7.82E+22 

7116. 

HN-ME 

61.2 

0.34 

50. 

2.7 

5.9 

3.36 

3.14E+22 

2861. 

KN-UT 

95.3 

0.69 

28. 

2.0 

3.15 

1.81 

3.73E+22 

3394. 

LC-NM 

93.6 

0.59 

53. 

2.2 

3.0 

1.6 

7.28E+22 

6628. 

OO-NW 

36.7 

0.46 

70. 

2.6 

5.8 

3.3 

3.86E+22 

3512. 

RK-ON 

76.6 

0.37 

56. 

2.7 

5.64 

3.47 

4.37E+22 

3982. 

SB-GR 

25.6 

0.55 

252.4 

2.6 

5.8 

3.3 

5.87E+22 

5343. 

SAPHIR 

BL-WV 

72.6 

0.37 

178. 

2.7 

6.0 

3.46 

1.32E+23 

12031. 

BR-PA 

70.2 

0.53 

206. 

2.7 

6.0 

3.46 

1.55E+23 

14105. 

HN-ME 

61.2 

0.34 

77. 

2.7 

5.9 

3.36 

4.85E+22 

4420. 

KN-UT 

95.3 

0.69 

52. 

2.0 

3.15 

1.81 

7.00E+22 

6374. 

LC-NM 

93.6 

0.59 

17o. 

2.2 

3.0 

1.6 

2.40E+2? 

21870. 

RT-NM 

89.8 

0.70 

334. 

2.6 

5.8 

3.3 

6.76E+23 

61521. 

RY-ND 

82.4 

0.57 

237. 

2.6 

5.8 

3.3 

1.94E+23 

17694. 

VO-IO 

78.9 

0.49 

196. 

2.6 

5.8 

3.3 

1.55E+23 

14103. 

WO-AZ 

95.1 

0.80 

484, 

2.6 

5.8 

3.3 

9.91E+23 

90171. 

M0  =  4itpa2  y,  was  assumed  to  relate  the  static  moment,  M0,  to  the  RDP  level,  4/„. 
Several  statistics,  shown  in  Table  4,  are  provided  for  the  "network"  estimate  of  M0,  including  a 


straight  average,  the  rms  value,  and  the  geometric  mean.  Following  McLaughlin  (1986)  we 
consider  the  rms  value  to  be  the  best  estimate  of  the  absolute  static  moments.  However,  the 


ratio  of  the  geometric  mean  (logarithmic  average)  should  be  a  better  estimate  of  the  relative 
source  size  of  the  two  events.  The  '■tations  RT-NM,  and  VV0-AZ  have  large  amplitudes  that 
may  be  due  to  deep  soft  sediment  amplification.  The  "network"  M0  estimates  appear  to  be  log- 
normally  distributed  and  the  variance  dominated  by  focusing-defocusing  amplitude  fluctuations 
with  a  C  for  !og10(Mo)  of  about  0.24  magnitude  units  for  RL'BIS  and  0.39  for  SAPHIR.  The 
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four  sites  that  recorded  both  RUBIS  and  SAPHIR  give  source  size  ratios  of  0.29,  0.65,  0.54, 


and  0.31  for  BL-WV,  HN-ME,  KN-UT,  and  LC-NM,  respectively.  The  average  log-ratio  of 


these  four  stations  is  -0.37.  This  is  consistent  with  the  log-ratio  of  the  yields  of  the  two 


events,  -0.36. 


TABLE  4. 


_ LRSM  SPECTRAL  MOMENT  ESTIMATES  1024  dyne-cm 

EVENT  MEAN  RMS  GEOMETRIC  MEAN  MEDIAN 
RUBIS  0.064(0.01)  0.072  0.056  0.059 

SAPHIR  0.30(0.10)  0.42  0.19  0.16 


Der  et  al.  (1982)  have  proposed  frequency  dependent  attenuation  models,  t*(f),  for  tec¬ 


tonic  and  shield  regions.  Since  t*  *»  t*  +  f-77-,  the  spectral  correction  made  for  frequency 

df 


?  dt* 

independent  t*  may  be  in  error  by  the  factor  exp(-Jt  r  — — ).  The  Der  et  al.  (1982)  tectonic 

df 


model  for  t*(f)  has  a  slope  at  1  Hz  of  about  -0.2  sec/Hz.  Therefore,  the  estimates  of  M0  and 


could  be  larger  by  a  factor  of  about  1.9.  This  factor  is  largely  uncertain  due  to  the  uncer¬ 


tainty  in  the  frequency  dependence  of  Q(f)  in  the  mantle  on  a  regional  basis  in  the  0.5  to  2.0 


Hz  bandwidth. 


DECONVOLUTIONS  AT  EKA  AND  YKA 


The  method  of  Shumway  and  Der  (1985)  was  used  to  deconvolve  the  array  data  for 


EMERAUDE,  GRENAT,  RUBIS,  and  SAPHIR  at  EKA  and  YKA.  A  constant  t*  -  0.45 


seconds  was  assumed  (Der  et  al.,  1985a  and  1985b).  Estimates  for  the  equivalent  far-field 


seismic  source  displacement  time  functions  are  shown  in  Figure  4  for  EKA  and  YKA.  The 


estimates  for  SAPHIR  differ  between  EKA  and  YKA  with  the  presence  of  a  negative  pulse 


following  the  positive  pulse  at  EKA  and  the  absence  of  such  a  pulse  at  YKA.  RUBIS  was 
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poorly  recorded  at  YKA.  The  absence  of  the  negative  pulse  (pP?)  following  the  main  SAPHIR 
P  pulse  at  YKA  as  well  as  the  variation  between  EKA  and  YKA  for  all  of  the  events  suggest 
that  the  topography  in  the  immediate  vicinity  of  the  explosions  may  have  influenced  the  later 
arrivals.  The  azimuths  to  EKA  and  YKA  from  SAPHIR  and  RUBIS  are  shown  in  Figure  2. 


Estimates  of  the  integrated  area  of  the  positive  pulses  are  shown  in  Figures  5A,  B,  C  and 
D  and  tabulated  in  Table  5.  These  measurements  are  similar  to  the  measurements  of  pulse  area 
made  from  deconvolved  records  used  by  Douglas  et  al.  (1986)  to  estimate  source  size  for  the 
three  Amchitka  explosions.  Two  methods  were  used.  In  the  first  method,  area  was  estimated 
within  the  pulse  defined  by  the  boxes  indicated  in  Figure  5A  and  B.  In  the  second  method,  the 
displacement  far-field  source  time  function  estimate  was  integrated  and  the  causal  step  was 
estimated  as  indicated  in  Figure  5C  and  5D. 

Also,  the  spectral  intercept  was  estimated  from  the  spectra  at  EKA  as  with  the  LRSM 
data.  The  spectral  estimate  for  was  identical  to  the  average  spectral  estimate  for  made 
at  EKA  for  RUBIS.  The  M0  estimates  for  RUBIS  and  SAPHIR  fall  within  the  population  of 
values  from  the  LRSM  network  of  stations,  Table  3.  The  network  scatter  far  exceeds  any  for¬ 
mal  error  in  the  estimates  made  at  any  given  station  within  the  network. 


TABLE  5. 


M0  1024 

dyne-cm 

nr  ! 

EVENT 

EKA 

YKA 

| 

SAPHIR 

~0705~ 

0.03  ' 

L  4600 

“2700  ; 

RUBIS 

0.03 

!  2700 

— 

GRENAT 

0.008 

0.004 

730 

360 

EMERAUDF. 

0.01 

0  00? 

90n 

730 

LRSM  TIME  DOMAIN  DECONVOLUTIONS 
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The  waveform  data  for  RUBIS  and  SAPHIR  were  deconvolved  to  ground  displacement  at 
several  of  the  LRSM  stations  by  removing  the  effects  of  the  instrument  and  an  estimated 
attenuation  operator.  The  Azimi  et  al.  (1968)  constant  t*  operator  was  used.  In  order  to  limit 
the  bandwidths  of  the  deconvolution,  a  filter  was  applied  to  each  deconvolution  that  maximizes 
the  useful  signal  bandwidth.  The  time  domain  response  of  this  filter  is  referred  to  as  the  reso¬ 
lution  kernel.  The  filter  was  of  the  form 

F(f)  =  - 1 -  - - - 

_0__  l+exp(-27C(f-fc)t*) 

|H(f)i2 

where  H(f)  is  the  complex  response  of  the  LRSM  instrument,  0  is  a  chosen  signal -to-noise 
level  and  fc  is  a  high  frequency  cut-off.  The  filter  is  nearly  flat  for  frequencies  below  fc  and 
for  the  bandwidth  of  frequencies  for  which  |  H(f)  | 2  >  0.  The  filter  serves  to  limit  the  decon¬ 
volution  of  the  high  and  low  frequencies  where  signal-to-noise  ratio  is  poor.  Typical  values  of 
fc  were  4.0  Hz  for  low  ?  stations  like  RK-ON  and  3.0  Hz  for  high  t*  stations  such  as  WO- 
AZ  or  KN-UT. 

As  with  the  EKA  and  YKA  data,  the  P  wave  signal  moment  was  estimated  from  the  time 
domain  deconvolved  pulses  and  used  to  estimate  the  explosion  moment,  M0.  The  results  in 
Table  6  are  in  good  agreement  with  the  spectral  estimates  of  moment  presented  in  Table  3. 
The  summary  statistics  are  given  in  Table  7  for  the  LRSM  temporal  moment  estimates  together 
with  the  EKA  and  YKA  estimates.  The  principal  difference  between  these  LRSM  deconvolu¬ 
tions  and  the  deconvolutions  of  Shumway  and  Der  (1985)  which  we  applied  to  the  EKA  and 
YKA  data  is  that  the  multichannel  deconvolution  estimates  the  local  station  site  responses. 
The  LRSM  deconvolutions  have  not  had  estimates  of  the  local  station  site  response  removed. 
The  deconvolutions  shown  in  Figures  6A  and  6B  are  bandlimited  estimates  of  the  explosion 
source  convolved  with  any  propagation  complexities  and  therefore  only  those  characteristics 
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that  are  similar  across  the  network  may  be  reliably  attributed  to  the  source. 


TABLE  6. 


RUB1S  AND  SAPHIR  LRSM  TEMPORAL  MOMENT  ESTIMATES 

STATION 

A 

t* 

P 

a 

P 

Mo 

A  oo 

deg 

S 

nm-s 

gm/cc 

km/s 

km/s 

dyne -cm 

m 

RUB  IS 

BL-WV 

72.6 

0.35 

71. 

2  7 

6.0 

3.46 

5.26E+22 

4700 

DH-NY 

66.7 

0.45 

88. 

2.7 

6.0 

3.46 

6.34E+22 

5772 

EB-MT 

78.1 

0.65 

204. 

2.6 

5.8 

3.3 

1 .61  Ei-23 

14643 

EU-AL 

79.4 

0.45 

185. 

1.9 

2.7 

1.4 

8.74E+22 

7950 

FR-MA 

86.3 

0.55 

160. 

2.6 

5.8 

3.3 

1.26E+23 

12346 

GI-MA 

84.4 

0.45 

88. 

2.6 

5.8 

3.3 

7.34E+22 

6681 

HL-ID 

92.2 

0.50 

66. 

2.5 

5.2 

2.96 

1.13E+23 

10300 

HN-ME 

61.2 

0.50 

92. 

2.7 

5.9 

3.36 

5.79E+22 

5275 

KN-UT 

95.3 

0.70 

30. 

2.0 

3.15 

1.81 

4.01E+22 

3649 

LC-NM 

93.6 

0.60 

63. 

2.2 

3.0 

1.6 

8.59E+22 

7819 

RK-ON 

76.6 

0.35 

92 

2.7 

5.64 

3.47 

7.17E+22 

6517 

SAPHIR 

BL-WV 

72.6 

0.35 

202. 

2.7 

6.0 

3.46 

1.50E+23 

13653 

BR-PA 

70.2 

0.50 

225. 

2.7 

6.0 

3.46 

1.69E+23 

15391 

HN-ME 

61.2 

0.50 

240. 

2.7 

5.9 

3.36 

1.51E+23 

13761 

HN-ME 

61.2 

0.50 

(511.) 

2.7 

5.9 

3.36 

(3.23E-I-23) 

(29300) 

KN-UT 

95.3 

0.70 

47. 

2.0 

3.15 

1.81 

6.28E+22 

5717 

LC-NM 

93.6 

0.60 

170. 

2.2 

3.0 

1.6 

2.32E+23 

21100 

RT-NM 

89.8 

0.70 

286. 

2.6 

5.8 

3.3 

5.78E+23 

52554 

RY-ND 

82.4 

0.60 

334. 

2.6 

5.8 

3.3 

2.74E-23 

24936 

VO-IO 

78.9 

0.50 

241.2 

2.6 

5.8 

3.3 

1.90E-23 

17324 

WO-AZ 

95.1 

0.80 

412. 

2.6 

5.8 

3.3 

8.44E-23 

76757 

TABLE  7. 

_ TEMPORAL  MOMENT  ESTIVATES  10W  dync-cm 

EVENT  MEAN  RMS  ciEOMETRic'MEAhT MEDIAN 

RUB  IS  0.08 1(0. 01)  6.089  0.07T"'  " "  0.07  3" 

SAPHIR  0.25(0.07)  0.34  0.16  0.17 


The  temporal  estimates  are  statistically  better  (smaller  dispersion)  than  the  spectral  esti¬ 
mates  and  the  ratio  of  the  si7e  of  the  two  events  Is  closer  to  the  ratio  of  the  yields  as  given  by 
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Marshall  et  al.  (1979).  The  average  log-ratio  of  RUBIS  to  SAPHIR  at  the  five  sites,  BL-WV, 
HN-ME,  KN-UT,  LC-NM,  and  EKA  is  -0.34.  This  is  in  excellent  agreement  with  the  log-ratio 
of  the  yields.  If  we  assume  that  the  moment  of  SAPHIR  is  2.3  times  the  moment  of  RUBIS 
then  we  arrive  at  the  summary  statistics  of  Table  8. 


TABLE  8. 


TEMPORAL  MOMENT  ESTIMATES  1024  dyne-cm 
-CON  STRAINED- 


EVENT 

MEAN 

RMS 

GEOMETRIC  MEAN 

MEDIAN 

RUBIS 

0.094(0.02) 

0.12 

0.072 

0.073 

SAPHIR 

0.22(0.04) 

0.28 

0.17 

0.17 

The  deconvolved  waveforms  are  shown  in  Figures  6A  and  6B.  The  waveforms  are 
shown  in  order  of  decreasing  azimuth  from  the  source.  Several  interesting  features  in  the  later 
arrivals  can  be  correlated  between  stations  and  are  probably  due  to  near  source  scattering.  The 
RUBIS  records  at  HL-ID  and  FR-MA  show  a  positive  late  pulse  (labeled  "A")  that  is  remark¬ 
ably  similar  between  the  two  stations.  A  similar  late  secondary  phase  can  be  seen  on  other  sta¬ 
tion  records  (labeled  "(A)”).  The  RUBIS  records  at  KN-UT,  HN-ME,  and  LC-NM  are  also 
similar  in  appearance  (azimuths  of  315,  311,  and  308  degrees,  respectively).  These  similarities 
are  probably  not  due  to  common  station  effects  and  therefore,  may  be  attributed  to  near  source 
complexity. 

SAPHIR  at  HN-ME  shows  a  broad  P  pulse  and  therefore  a  large  time  domain  moment 
estimate  in  contrast  to  the  spectral  estimate  (in  parentheses  in  TABLE  6).  A  better  estimate  at 
HN-ME  may  be  to  ignore  this  broadening  since  the  initial  P  wave  is  simple  at  all  the  other  sta¬ 
tions.  This  may  be  due  to  a  secondary  arrival  prolonging  the  SAPHIR  P-wave  pulse  at  HN- 
ME.  Since  this  is  isolated  to  HN-ME  it  may  be  a  station  effect,  although  the  RUBIS  record 
may  not  show  the  same  feature. 
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DISCUSSION 

It  is  of  interest  to  compare  these  teleseismic  RDP  values  with  free-field  measurements  of 
RDP  for  the  three  US  granite  explosions,  PILEDRIVER,  SHOAL,  and  HARDHAT.  Werth 
and  Herhst  (1963)  present  data  for  HARDHAT  (5KT,p-2.69gnVcc,a=*4  8km s,h«286m)  with 
peak  values  of  RDP  near  4700  m3  and  a  static  level  of  about  2500  m\  Murphy  (1978) 
presents  data  for  PILEDRIVER  (62KT,p=2.67gm/cc,ot=5.8km,s,P=3.45km/s,h=457m)  with 
peak  levels  of  30000  to  50000  m3  and  static  levels  of  about  20000  m\  for  SHOAL 
(13KT,p=2.55gm/cc,a=5.5km/s,P=3.0km/s,h=367m)  with  peak  RDP  levels  of  7000  m3  and 
static  levels  of  about  4000  m\  and  for  HARDHAT 
(5.9KT,p=2.67gnVcc,a=5.5km/s,p=3.25km/s,h-290m)  with  peak  RDP  levels  of  3000  to  5000 
m3  and  static  levels  of  about  2300  m3. 

The  static  PILEDRIVER  moment  of  0.22  lO34  dyne-cm  (Murphy,  1978)  when  scaled  to 
the  RUBIS  yield  is  0.18  1024  dyne-cm.  This  is  1.5  times  the  RUBIS  rms  moment  estimate  of 
Table  8.  The  discrepancy  is  significantly  larger  if  the  mean  or  geometric  mean  of  the  RUBIS 
moment  is  used  for  comparison,  or  if  the  peak  PILEDRIVER  RDP  levels  are  used  for  com¬ 
parison.  Comparison  of  the  teleseismic  moment  estimates  with  the  scaled  free-field  RDP  levels 
from  SHOAL  or  HARDHAT  give  much  the  same  result.  The  French  Sahara  teleseismic  rms 
moment  estimates  are  between  30  and  50%  too  low  when  compared  to  the  measured  free-field 
RDP  levels  of  US  granite  explosions.  The  measured  free-field  granite  RDP  levels  are  roughly 
consistent  with  each  other,  while  Heuze  (1983)  and  Schock  rt  al  (1972)  tind  no  significant 
difference  between  the  granites  of  Climax  Stock,  NTS.  and  1  aourirt  Tam  Afe’la  Massif, 
Algeria.  However,  the  free-field  RDP  levels  are  clearlv  inconsistent  with  the  teleseismic  obser¬ 
vations  of  the  equivalent  moments  for  RUBIS  and  SAPH1R  unless  additional  attenuation  is 
introduced  into  either  the  near-field  or  along  the  teleseismic  path 
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Several  hypotheses  may  be  presented  to  explain  the  discrepancy  between  the  two  data 
sets.  First  of  all,  the  free-field  RDP  data  were  all  collected  within  the  non-linear  zone  of 
ground  motion  and  may  over  estimate  the  far-field  equivalent  explosion  moment.  Only 
attenu  ition  measurements  at  high  strain  rates  in  granite,  can  directly  test  this  hypothesis.  Some 
additional  broadband  attenuation  is  introduced  by  reflections  and  conversions  at  crustal  and 
mantle  interfaces,  but  it  is  unlikely  that  this  amounts  to  between  30  and  50%  broadband  reduc¬ 
tion  in  amplitude.  Similarly,  McLaughlin  and  Anderson  (1984)  suggest  that  crustal  scattering 
may  contribute  to  pulse  attenuation  but  may  not  be  detectable  by  spectral  t*  measurements. 
However,  their  evidence  suggests  that  the  mechanism  is  not  significant  at  frequencies  below  1 
dt* 

Hz.  Finally,  if  f  — —  =  -0.2  sec  in  the  0.5  to  2  Hz  bandwidth,  then  the  teleseismic  moment 
df 

estimates  presented  are  too  low  and  a  factor  of  2  correction  would  be  consistent  with  the  static 
free-field  granite  RDP  measurements  made  for  PILEDRIVER,  SHOAL,  and  HARDHAT. 
Several  such  frequency  dependent  t*(f)  models  have  been  proposed  (Der  et  al.,  1982;  Bache  et 
al.,  1985;  Taylor  et  al.,  1986;  Choy  and  Cormier,  1986;  and  Der  et  al.,  1986). 

CONCLUSIONS 


We  have  estimated  the  far-field  equivalent  seismic  sources  for  four  French  Sahara  explo¬ 
sions  contained  in  granite  from  short  period  seismic  recordings,  using  two  time  domain  decon¬ 
volution  methods  and  one  spectral  estimation  method.  Uncertainties  in  source  size  are  largely 
due  to  the  uncertainties  in  the  frequency  dependence  of  Q  in  the  mantle,  and  the  amplitude 
fluctuations  due  to  focusing-defocusing  along  the  teleseismic  path.  The  uncertainties  in  t*(f) 
near  1  Hz  are  systematic  errors  while  the  focusing-defocusing  errors  are  more  likely  random  in 

dt* 

nature.  The  most  serious  of  these  two  errors  is  probably  the  uncertainty  in  the  value  of  — 
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dt* 

near  1  Hz.  A  value  of -  of  about  -0.2  sec/Hz  at  1  Hz  would  roughly  double  the  moment 

df 

estimates.  The  focusing-defocusing  introduces  random  error  into  the  netwom  but  may  be  over¬ 
come  by  the  averaging  over  many  stations.  Focusing-defocusing  patterns  can  be  discerned  and 
corrected  from  mb  measurements  from  WWSSN  stations  or  the  ISC  reported  magnitudes. 

P  wave  teleseismic  t*  estimates  for  the  test  site  are  in  agreement  with  classification  of  the 
Ahaggar  region  as  a  "hot  spot".  The  values  of  t*  estimated  from  LRSM,  EKA,  and  YKA  data 
are  in  agreement  with  the  Der  et  al.  (1985b)  estimate  of  0.35  sec  but  favor  a  slightly  lower 
value  near  0.31  sec.  The  Der  et  al.  (1985b)  t*  estimate  for  NTS  is  0.34  sec,  therefore  there  is 
no  evidence  for  significant  attenuation  bias  between  Ahaggar  and  NTS. 

The  deconvolution  technique  of  Shumway  and  Der  (1985)  has  been  shown  to  yield 
explosion  size,  M0  or  4,ao  values,  consistent  with  spectral  estimates.  SAPHIR  recorded  at 
YKA  may  show  evidence  of  destructive  interference  of  pP  due  to  topographic  effects. 

Deconvolutions  of  the  LRSM  records  reveal  similar  variability  of  the  radiated  P 
waveforms  from  SAPHIR  and  RUBIS,  Particularly  interesting  are  the  late  arrivals  correlated 
between  stations  some  distance  apart  but  of  similar  takeoff  angle.  The  topography  at  the 
French  southern  Sahara  test  site  could  be  responsible  for  these  large  variations. 

Comparison  of  the  various  methods  used  in  this  work  reveals  that  the  spectral  and  decon¬ 
volution  methods  make  consistent  estimates  of  explosion  size.  While  the  spectral  slope  and 
intercept  method  may  be  used  to  estimate  the  t*  and  U0  for  the  pulse,  the  method  relies  on  a 
rough  estimate  of  the  source  spectral  shape.  The  deconvolution  methods  do  not  assume  any¬ 
thing  about  the  source  time  function  but  require  an  estimate  of  the  attenuation  parameters. 
Therefore  it  is  not  surprising  that  the  deconvolution  methods  that  use  the  spectral  t*  estimates 
give  results  consistent  with  the  low  frequency  spectral  intercept  method  However,  the  decon¬ 
volutions  reveal  a  great  deal  of  complexity  in  the  early  P-wave  coda  while  the  initial  P  wave 
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remains  very  simple.  Spectral  methods  average  the  P  waveform  and  may  be  susceptible  to  the 
interference  of  positive  and  negative  polarity  pulses  to  decrease  the  low  frequency  level  of  the 
P  waveform.  The  presence  of  these  interfering  negative  and  positive  pulses  can  then  be  sorted 
out  using  the  deconvolution  methods.  By  this  token  the  time  domain  moment  estimates  are  the 
better  estimates  and  justify  the  extra  computation  and  care  that  must  be  made  in  deconvolution 
of  bandlimited  signals. 
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FIGURE  CAPTIONS 

FIGURE  1.  Locations  of  LRSM  stations,  the  EKA,  and  YKA  arrays  centered  on  the  southern 
Algerian  test  site.  The  two  arrays,  EKA  (Eskdalemuir,  Scodand,  U.K.)  and  YKA  (Yellow¬ 
knife,  Northwest  Territories,  Canada)  are  indicated  with  solid  symbols. 

FIGURE  2.  Contour  map  of  Taourirt  Tan  Afella  from  Faure  (1972)  and  locations  for  several 
events  from  Duclaux  and  Michaud  (1970)  are  indicated.  The  contours  are  at  100  meter  inter¬ 
vals  and  the  dotted  line  shows  the  granite  outcrop.  The  elevation  of  the  plateau  around  the 
mountain  is  about  1100  meters.  Locations  for  SAPHIR  (S),  RUBIS  (R),  EMERAUDE  (Em), 
and  GRENAT  (G)  are  indicated.  The  azimuths  to  several  stations  are  shown  above  with 
labeled  arrows. 

FIGURE  3.  Spectra  corrected  for  instrument  and  RDP.  Slopes  between  0.5  and  3.0  Hz  are 
estimates  of  t*,  intercepts  are  estimates  of  D0.  Only  spectral  levels  with  an  estimated  signal- 
to-noise  power  ratio  greater  than  2  were  used. 

FIGURE  4.  Deconvolved  source  time  functions  at  the  EKA  and  YKA  arrays  for  SAPHIR, 
RUBIS,  GRENAT  and  EMERAUDE.  RUBIS  was  poorly  recorded  across  YKA. 

FIGURE  5A.  Equivalent  seismic  sources,  M(t)  for  SAPHIR,  RUBIS,  GRENAT  and 

EMERAUDE  at  the  EKA  array.  The  explosion  moment  inferred  from  the  P-pulse  area  is  indi¬ 
cated. 

FIGURE  5B.  Equivalent  seismic  sources,  M(t),  for  SAPHIR,  RUBIS,  GRENAT  and 

EMERAUDE  at  the  YKA  array.  The  explosion  moment  inferred  from  the  P-pulse  area  is  indi¬ 
cated. 

FIGURE  5C.  Equivalent  seismic  sources,  M(t),  for  SAPHIR,  RUBIS,  GRENAT  and 

EMERAUDE  at  the  EKA  array.  The  explosion  moment  inferred  from  the  P-pulse  area  is  indi¬ 
cated. 

FIGURE  5D.  Equivalent  seismic  sources,  M(t),  for  SAPHIR,  RUBIS,  GRENAT  and 

EMERAUDE  at  the  YKA  array.  The  explosion  moment  inferred  from  the  P-pulse  area  is  indi¬ 
cated. 

FIGURE  6A.  Deconvolved  waveforms  for  RUBIS  at  LRSM  stations.  The  causal  P-wave  pulse 
area  is  indicated  for  each  trace.  Traces  are  ordered  in  decreasing  azimuth  from  the  top.  Late 
positive  arrivals  that  may  be  correlated  between  stations  are  indicated  by  labels  "A"  and  "(A)". 

FIGURE  6B.  Deconvolved  waveforms  for  SAPHIR  at  LRSM  stations.  The  causal  P-wave 
pulse  area  is  indicated  for  each  trace.  Traces  are  ordered  in  decreasing  azimuth  from  the  top. 
The  HN-ME  waveform  shows  considerable  broadening  that  may  be  due  to  a  positive  secondary 
arrival  unique  to  the  HN-ME  station. 
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SECTION  3 

SCATTERING  FROM  NEAR-SOURCE  TOPOGRAPHY: 

TELESEISMIC  OBSERVATIONS  AND 
NUMERICAL  2-D  EXPLOSIVE  LINE  SOURCE  SIMULATIONS 

K.  L.  McLaughlin,  and  R.-S.  Jih 
Teledyne  Geotech  Alexandria  Laboratories 
314  Montgomery  Street 
Alexandria,  VA  22314 

INTRODUCTION 

Jih  et  al  (1986)  propose  methods  for  introduction  of  free  surface  topography  in  hetero¬ 
geneous  elastodynamic  two  dimensional  (2-D)  linear  finite  difference  (LFD)  calculations.  We 
have  applied  the.ve  methods  to  explosive  line  sources  for  simulation  of  far-field  "teleseismic"  P 
waves.  The  motivation  for  these  numerical  experiments  is  to  examine  the  extent  of  scattering 
from  near-source  topography  on  short-period  seismograms  used  for  explosion  yield  estimation. 
Several  important  test  sites  have  substantial  topography  within  a  few  (1  Hz  P  wave) 
wavelengths  of  the  source.  These  test  sites  include  the  French  southern  Sahara  test  site,  the 
Novaya  Zemlya  test  site  (Greenfield,  1971),  and  the  Degelen,  Eastern  Kazakh  test  site  (Rodean, 
1979). 

To  explore  the  variability  that  may  be  introduced  by  near-source  topography,  several  sim¬ 
ple  test  cases  have  been  investigated.  These  test  cases  •  e  of 

the  far-field  teleseismic  P  the 

;r>v  :  •  e. 
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As  a  worst  case,  the  southern  Algeria  test  site  located  within  Taourirt  Tan  Afella  Massif 
will  serve  as  a  model  for  some  more  specific  calculations.  This  granite  mountain  rises  from 
the  (1100  meter  elevation)  Ahaggar  plateau  to  a  peak  of  2100  meters  in  less  than  1  km.  The 
mountain  is  only  16  km2  and  has  slopes  that  locally  exceed  45°.  Several  tests  were  conducted 
at  this  site  in  southern  Algeria  from  1962  to  1966.  A  topographic  map  from  Duclaux  and 
Michaud  (1970)  is  shown  in  Figure  1,  with  locations  from  Faure  (1972)  (see  Table  1).  Con¬ 
tours  are  at  100  meter  intervals  and  the  granite  outcrop  is  indicated  by  the  dashed  line.  In 
recent  years  Blandford  and  Shumway  (1982),  Blandford  et  al  (1984)  and  McLaughlin  et  al 
(1986a)  have  conducted  investigations  of  the  mb  observed  from  these  explosions  at  WWSSN 
stations.  Furthermore,  available  Long  Range  Seismic  Measurements  (LRSM)  data  and  the  Esk- 
dalemuir,  Scotland  (EKA)  and  Yellow  Knife,  Northwest  Territories  (YKA)  arrays  have  been 
analyzed  (see  McLaughlin  et  al  1986b)  for  estimates  of  the  explosion  spectra,  moments,  and 
source  time  functions. 


TABLE  1. 

EVENT  DATA  FROM  FAURE  (1972) 

DATE 

ORIGIN  TIME 

LAT  (N) 

LONG  (E) 

EVENT 

1963/03/18 

10:02:00.351 

24.0414 

5.0522 

EMERAUDE 

1963/10/20 

13:00:00.011 

24.0355 

5.0386 

RUBIS 

1965/02/27 

11:30:00.039 

24.0587 

5.0312 

SAPHIR 

1965/12/01 

10:30:00.088 

24.0437 

5.0469 

TOURMALINE 

1966/02/16 

11:00:00.035 

24.0442 

5.0412 

G REN AT 

AN  EXAMPLE  TOPOGRAPHIC  CALCULATION 

Figure  2 A  shows  a  simple  topographic  profile.  The  elastic  medium  has  a  P-wave  velo¬ 
city  of  5  km/s.  and  the  line  source  locations  are  indicated  by  triangles  in  the  Figure.  The  cal¬ 
culations  are  performed  by  use  of  reciprocity.  A  broadband  plane  wave  (source  at  infinity)  is 
incident  upon  die  model  and  the  dilatation  time  history  is  measured  at  a  location  within  the 
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grid.  By  reciprocity,  this  time  history  is  the  far-field  broadband  displacement  pulse  at  infinity 
for  a  broadband  dilatational  source  within  the  model.  We  may  simulate  any  far-field  seismo¬ 
gram  by  the  removal  of  the  known  bandhmited  incident  P-wave  pulse,  and  convolution  of  a 
teleseismic  transfer  function  that  includes  attenuation,  receiver  crustal  structure,  instrument 
response,  and  model  source  time  function  (c.f.  Robinson,  1983).  A  single  run  with  an 
incident  planar  P  wave  then  yields  the  teleseismic  response  for  a  suite  of  source  locations 
within  the  model.  A  run  must  be  performed  for  each  takeoff  angle  (teleseismic  slowness) 
desired  by  introduction  of  a  planar  P  wave  with  the  desired  incidence  angle. 

Figure  3  shows  teleseismic  synthetic  seismograms  for  the  the  topographic  profile  in  Fig¬ 
ure  2A  for  a  takeoff  angle  of  0°.  The  synthetic  has  been  convolved  with  a  von  Seggern  and 
Blandford  (1972)  granite  explosion  source  time  function  (100  KT),  an  attenuation  operator  (t* 
=  0.45  sec),  and  an  LRSM  instrument  response.  All  of  the  seismograms  are  plotted  at  the 
same  scale  and  show  die  variability  of  the  P-waveforms  due  to  the  coastructive  and  destructive 
interference  of  the  linear  reflection  from  the  surface.  The  "be”  phase  of  the  P  wave  shows  the 
greatest  variability  due  to  the  elastic  P+pP  interference.  The  log-amplitudes  of  the  P  wave  "a”, 
"ab",  and  "max"  phases  in  Figure  3  are  plotted  in  Figure  2B.  The  iog(a)  variation  is  slight  as 
expected  while  the  log(ab)  and  log(max)  variations  reflect  considerable  variation  in  the  linear 
elastic  P+pP  interference. 


TAOURIRT  TAN  A  FELL  A  MASSIF 

A  north-to-south  cross  section  of  Taourirt  Tan  Afella  Massif  (P  wave  velocity  of  5.2 
km/s)  is  shown  in  Figure  4A  with  a  1.5  incident  pianat  P  wave.  A  planar,  band  limited  P 
wave  is  incident  upon  the  medium  from  below  and  the  dilatational  strain  is  measured  at  the 
source  locations.  The  far-field  di%p';vement  P  waveform  due  to  m  explosion  line  source  at 
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these  locations  is  then  inferred  from  this  transfer  function  by  reciprocity.  Figures  4B  through 
4E  illustrate  the  dilatational  field  at  intervals  of  0.2  second  as  the  P  wave  reflects  from  the  free 
surface.  Note  that  at  time  interval  0.4  second,  the  strong  constructive  and  destructive  interfer¬ 
ence  that  results  from  the  reflections  directly  beneath  the  peaks  of  the  ridge  (2-D  mountain). 
The  linear  pP  contribution  from  this  location  would  be  large  and  variable  for  a  receivers  at  a 
takeoff  angles  near  15°.  At  time  intervals  0.6  and  0.8  seconds  we  see  what  appear  to  be  Ray¬ 
leigh  waves  trapped  at  the  surface  as  well  as  dilatational  coda  behind  the  specularly  reflected  P 
wave. 

Source  locations  indicated  by  triangles  in  Figure  4A  were  used  to  infer  far-field  P 
waveforms  from  2-D  explosive  line  sources  shown  in  Figure  5A  and  5B.  The  waveforms  have 
been  convolved  with  von  Seggem  and  Blandford  (1972)  explosive  source  time  functions  (50 
and  120  KT),  an  attenuation  operator  (t*  -  0.45  sec),  and  an  LRSM  instrument  response 
Waveforms  are  shown  for  takeoff  angles  of  0,  5,  10,  15,  and  20  degrees.  These  two  source 
locations  are  labeled  "RUBIS"  and  "SAPHIR"  analogous  to  the  two  locations  (R  and  S)  of 
Figure  1.  Several  interesting  features  are  evident  in  these  2-D  numerical  simulations. 

For  the  location  "SAPHIR",  the  maximum  amplitude  decreases  with  increasing  takeoff 
angle  (see  Figure  5A).  The  elastic  pP  contribution  increases  in  such  as  way  as  to  reduce  the 
"b-c"  amplitude  with  increasing  takeoff  angle.  The  coda  amplitude  increases  with  increasing 
takeoff  angle  as  the  free-surface  interaction  becomes  more  complicated  and  the  conversions  of 
P-SV  energy  increase  at  the  free  surface.  The  "SAPHIR'  location  is  roughly  coincident  with 
the  focal  point  of  the  topographic  profile  and  the  pP  reflection  may  be  amplified  at  this  source 
location  for  some  narrow  range  of  takeoff  angles 

The  "RUBIS"  location  also  shows  an  increase  of  coda  amplitude  with  increasing  takeuft 
angle  (See  Figure  5B).  For  the  shallowest  takeoff  angle  of  20",  there  is  a  semndurv  arns.il 
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larger  than  the  main  P  wave.  The  "RUBIS"  ”b-c"  amplitude  does  not  show  as  strong  a  varia¬ 
tion  with  takeoff  angle.  The  topography  directly  above  the  "RUBIS"  location  is  less  steep  and 
the  pP  free  surface  interaction  is  only  distorted.  However,  the  topography  to  the  north  of  the 
"RUBIS"  location  introduces  scattered  arrivals  responsible  for  the  strong  coda.  The  "RUBIS" 
location  is  above  the  "focal  point"  of  the  topographic  profile  and  the  pP  reflection  may  be  de- 
focused  at  this  source  location. 

Spectral  ratios  of  the  "RUBIS" -to-"SAPHIR"  2-D  synthetics  showed  that  as  the  takeoff 
angle  increases,  "RUBIS"  becomes  enriched  in  high  frequencies  with  respect  to  "SAPHIR". 
The  modulation  of  the  spectra  makes  this  comparison  difficult  if  the  bandwidth  is  limited  as  it 
is  in  real  data  to  frequencies  below  4.0  Hz.  However,  the  spectral  ratio  of  the  two  events  at  a 
common  station  may  give  clues  as  to  which  event  received  the  largest  contribution  of  scattered 
energy  in  the  coda  due  to  the  topography  above  the  source. 

YKA,  EKA,  KN-UT,  HN-ME,  BL-WV,  and  LC-NM  recorded  both  RUBIS  and 
SAPHIR.  The  spectral  ratios  of  RUBIS-to-SAPHIR  are  shown  in  Figure  6.  The  results  from 
the  spectral  ratios  indicate  that  the  "average"  slope  between  the  two  events  over  the  stations  is 
near  zero,  and  that  some  of  the  stations  show  positive  slopes  while  others  show  negative 
slopes.  The  bandwidths  are  in  most  cases  limited  to  between  0.5  and  3.0  Hz.  The  azimuths  of 
these  stations  span  about  45°  (due  north,  EKA,  to  northwest,  BL-WV,  see  Figure  1).  The  data 
indicate  that  RUBIS  and  SAPHIR  have  similar  average  frequency  content  and  that  the  effect  of 
scattering  may  vary  rapidly  with  azimuih  in  the  real  3-D  world. 

COMPARISONS  WITH  DECONVOLVED  DATA 

Four  events  EMERAUDE,  GRENAT,  RLB1S,  and  SAPHIR  were  deconvolved  at  the 
EKA  and  YKA  arrays  using  the  method  of  Shumway  and  Der  (1985).  The  results  are  shown 
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in  Figure  7.  The  equivalent  seismic  sources  for  the  events  are  shown  side-by-side  for  the  two 
arrays  and  a  comparison  can  be  made  for  EMERAUDE,  GRENAT,  and  SAPHIR  (RUBIS  was 
poorly  recorded  at  YKA).  The  initial  P  waveform  is  similar  for  each  event  at  the  two  arrays 
while  the  differences  are  greatest  beginning  about  0.5  sec  after  the  initial  P  wave.  SAPHIR 
shows  the  clearest  differences  between  the  two  arrays;  the  negative  pulse  following  the  positive 
pulse  is  strong  at  EKA  but  absent  at  YKA.  EMERAUDE  and  GRENAT  show  differences  in 
timing  of  the  negative  (and  small  positive)  pulses  following  the  initial  P  wave.  These 
differences  argue  in  favor  of  a  scattering  model  for  the  variations  in  the  free  surface  interaction 
for  these  events. 

Deconvolutions  of  available  LRSM  data  for  SAPHIR  and  RUBIS  are  presented  in  Fig¬ 
ures  8A,  and  8B.  The  waveforms  show  bandlimited  estimates  of  the  ground  displacement  at 
each  station  for  the  two  events.  The  instrument  response  and  a  constant  t*  attenuation  operator 
(Futterman,  1962;  Azimi  et  al.,  1968)  have  been  removed  from  the  digital  waveforms  and  a 
bandlimiting  filter  has  been  applied  to  maximize  the  broadband  signal-to-noise  ratio  of  the  data. 
The  data  is  generally  limited  to  frequencies  between  0.25  and  4.0  Hz.  Lower  Q  stations  are 
limited  to  frequencies  below  3.0  Hz.  Details  of  the  deconvolution  procedure  are  presented  in 
McLaughlin  et  al  (1986b). 

The  displacement  waveforms  of  Figures  8A  and  8B  are  arranged  with  increasing  source- 
to-receiver  azimuth  from  top  to  bottom.  The  area  under  the  initial  causal  P  waveform  is 
shaded.  There  are  several  characteristics  of  the  waveforms  that  are  readily  apparent. 

First,  the  initial  displacement  P  wave  is  rather  simple  with  a  few  exceptions.  Because  of 
the  bandlimited  nature  of  the  deconvolutions,  features  with  time  constants  shorter  than  0.3 
seconds  are  probably  not  resolvable.  In  fact  the  variability  of  the  initial  P  waveform  with  few 
exceptions  could  be  attributed  to  variations  in  the  resolution  kernels  of  the  deconvolutions. 
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Second,  there  are  several  positive  pulses  (labeled  "A"  and  "(A)")  that  can  be  correlated 
from  station-to-station  with  similar  takeoff  angles  and  azimuths.  For  example,  the  two  stations 
GI-MA  (Montana)  and  HL-ID  (Idaho)  show  a  similar  positive  arrival  labeled  "A”  from  the 
explosion  RUBIS.  Although  RK-ON  is  at  nearly  the  same  azimuth  as  these  two  stations,  RK- 
ON  is  at  a  greater  distance  on  the  focal  sphere  from  either  of  the  two  stations  because  of 
dependence  of  takeoff  angles  on  epicentral  distance.  If  these  arrivals  are  generated  near  the 
source,  they  would  indicate  significant  direction  properties  of  the  P  coda  generation  process. 

Third,  there  is  little  or  no  consistent  simple  "pP“  at  the  majority  of  stations  for  either 
SAPHIR  or  RUBIS.  Those  stations  that  do  show  negative  excursions  have  generally  longer 
period  and  possibly  multiple  negative  arrivals  following  the  first  positive  polarity  P  pulse  by 
0.5  to  1.0  seconds.  Given  the  high  P-wave  velocity  (about  5.8  km^s)  of  the  granite  mountain, 
Taourirt  Tan  Affela,  the  "linear  elastic”  pP-P  times  should  be  considerably  less  than  0.5 
seconds  for  both  explosions. 

Several  fundamental  differences  exist  between  2-D  and  3-D  wave  propagation  that 
prevent  quantitative  comparison  of  observations  and  synthetics.  The  most  significant  of  these 
differences  may  be  that  2-D  Rayleigh  waves  from  line  sources  do  not  have  geometrical  spread¬ 
ing  and  that  2-D  P-waves  attenuate  proportional  to  rather  than  — .  Therefore  we  should 

vr  r 

expect  that  the  teleseismic  P-coda  from  line  sources  in  a  2-D  structure  will  over  estimate  the 
effects  of  3-D  scattering.  Because  of  this  we  can  make  only  qualitative  comparisons  of  the 
observations  with  the  synthetic  calculations.  Given  this  caveat,  we  can  draw  some  parallels 
between  the  synthetics  and  the  deconvolved  LRSM  records. 

P  coda  elements  common  to  different  stations  with  similar  takeoff  angle  and  azimuth  may 
be  attributed  to  the  source  region.  The  longer  period  pulses  that  follow  the  RUBIS  P  wave  by 
about  1  second  are  seen  at  GI-MA.  KN-UT,  HN-MF.,  and  l.C-NM.  Similar  long-period 
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positive  polarity  secondary  pulses  may  be  seen  on  the  synthetics  of  Figures  9A  and  9B.  The 
synthetic  positive  polarity  secondary  arrivals  appear  to  be  scattered  Rayleigh-to-P  from  the 
topography  as  Rayleigh  waves  excited  on  the  free  surface  are  scattered  at  topographic  slope 
changes  and  contribute  to  P  coda. 

Another  characteristic  of  the  observed  data  is  the  variability  of  any  pP  arrivals.  The 
linear  elastic  pP  pulse  is  expected  to  follow  the  P  wave  much  earlier  than  observed.  No  such 
negative  polarity  pulse  is  consistently  observed  within  0.5  seconds  on  the  majority  of  the 
waveforms.  Some  displacement  waveforms  suggest  greatly  delayed  or  multipathed  pP  pulses 
(SAPHIR  at  BL-WV,  BR-PA,  HN-ME,  and  VO-IO,  and  RUBIS  at  EU-AL,  BL-WV,  DH-NY, 
HN-ME,  KN-UT,  FR-MA,  GI-MA,  RK-ON,  and  HL-ID).  This  same  characteristic  may  be 
seen  in  the  synthetics  at  various  azimuths  and  takeoff  angles.  The  2-D  linear  elastic  pP 
reflection  can  be  focused  (or  defocused)  by  constructive  (or  destructive)  reflections  from  the 
topography.  This  linear  elastic  P+pP  interference  introduces  fluctuations  in  the  "be"  swing  of 
the  teleseismic  P-wave  seismograms. 

Non-linear  spall  could  also  produce  these  same  two  characteristics  seen  in  the  data  and 
the  2-D  synthetics.  Non-linear  spall  could  produce  multiple,  lengthened  and  delayed  "pP" 
arrivals  as  well  as  other  secondary  phases.  However,  the  elastic  2-D  linear  finite  difference 
calculations  indicate  that  at  least  some  of  these  effects  may  be  due  to  steep  topographic  scatter¬ 
ing.  The  spall  mechanisms  may  not  be  required  to  explain  the  lack  of  clear  pP  arrivals  for  the 
RUBIS  and  SAPHIR  teleseismic  displacement  waveforms. 

WWSSN  MAGNITUDE  VARIATIONS 

In  order  to  address  the  question  as  to  whether  the  short-period  m*,  magnitude  distribution 
has  any  evidence  of  scattering,  we  have  examined  the  measurements  of  mb(Pa),  mb(Pb),  and 
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mb(Pmax)  made  at  WWSSN  stations  for  RUBIS,  SAPHIR,  GRENAT,  and  TOURMALINE. 
The  "Pa"  measurement  is  the  zero-to-peak  amplitude  of  the  first  motion  of  the  short  period 
seismogram,  "Pb"  is  1/2  of  the  a-to-b  peak-to-peak  measurement.  "Pmax"  is  1/2  the  maximum 
peak-to-peak  measurement  in  the  first  5  seconds  of  the  P  arrival. 

The  magnitudes  were  estimated  in  a  maximum  likelihood  sense  (Ringdal,  1976)  and  are 
listed  in  Table  2.  The  data  is  corrected  for  the  biasing  effects  of  non-detection  and/or  clipping. 
In  the  case  of  non-detection,  the  noise  level  is  measured  at  the  expected  arrival  time  for  the  P 
wave  and  used  as  the  detection  threshold.  In  the  case  of  clipping,  a  conservative  estimate  was 
made  of  the  maximu  visible  peak-to-peak  deflection  on  the  seismogram  and  used  as  the  clip¬ 
ping  level.  The  station  corrections  of  Ringdal  (1986)  and  the  distance  corrections  of  Veith  and 
Clawson  (1972)  were  applied  to  the  individual  station  magnitudes.  Uncertainties  in  the  event 
magnitude  were  estimated  by  use  of  a  bootstrap  procedure  (Efron,  1981).  PILEDRIVER  and 
SHOAL  are  included  for  comparison;  these  are  the  only  US  granite  explosions  for  which 
WWSSN  data  is  available. 


TABLE  2 

MAXIMUM  LIKELIHOOD  MAGNITUDE  ESTIMATES 

EVENT 

mb(Pa) 

mb(Pb) 

mb(Pmax) 

RUBIS 

4.76(.08) 

5. 14(.06) 

5.42(.06) 

SAPHIR 

5. 1 8(.06) 

5.49(.07) 

5.75(.06) 

GRENAT 

4.04(.  12) 

4.41(.07) 

4.69(.07) 

TOURMALINE 

3.99(.20) 

4.3 1  (.  10) 

4.60(.07) 

PILEDRIVER 

4.93(.05) 

5.21(.04) 

5.46(.05) 

SHOAL 

4.09(.  10) 

4.32(.08) 

4.63(.08) 

As  expected,  the  smallest  events  are  the  most  uncertain.  The  French  Sahara  magnitude 
distributions  have  standard  deviations  of  a  single  observation  of  0  40(0.05)  magnitude  units. 
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The  mb(Pa),  mb(Pb),  and  mb(Pmax)  are  correlated  across  the  network  such  that  the  loga¬ 
rithmic  ratios,  log(Pmax/Pa),  log(Pb/Pa),  and  log(Pmax/Pb)  are  well  defined.  These  ratios  are 
estimated  in  a  maximum  likelihood  sense  and  listed  in  Table  3.  The  standard  deviation  of  the 
log-ratio  distributions,  a,  are  also  listed.  The  estimates  of  log(Pmax/Pa)  and  log(Pb/Pa)  for 
GRENAT  and  TOURMALINE  may  not  be  as  well  determined  as  those  for  SAPHIR  and 
RUBIS. 


TABLE  3 

MAXIMUM  LIKELIHOOD  LOG-RATIO  ESTIMATES 

EVENT 

log(Pmax/Pa) 

log(Pb/Pa) 

log(Pmax/Pb) 

a 

o 

a 

RUBIS 

0.55 

0.27 

0.26 

0.18 

0.23 

0.17 

SAPHIR 

0.53 

0.28 

0.25 

0.19 

0.21 

0.19 

GRENAT 

0.67 

0.45 

0.36 

0.42 

0.24 

0.16 

TOURMALINE 

0.70 

0.45 

0.32 

0.34 

0.25 

0.17 

PILEDRIVER 

0.52 

0.28 

0.25 

0.19 

0.22 

0.18 

SHOAL 

0.49 

0.22 

0.31 

0.16 

0.31 

0.16 

Comparison  of  the  o  values  listed  in  Table  3  with  the  a  -  0.4  value  found  for  the 
mb(Pa),  mb(Pb),  and  mb(Pmax)  distributions  give  estimates  of  the  normalized  correlations  of 
log(Pa),  log(Pb),  and  log(Pmax).  log(Pmax)  and  log(Pb)  are  80%  correlated  across  the  net¬ 
work  for  all  of  the  explosions,  log(Pb)  and  log(Pa)  are  70%  correlated  across  the  network,  but 
log(Pa)  and  log(Pmax)  are  only  40%  correlated  across  the  network  for  RUBIS  and  SAPHIR. 
Since  the  log(Pa)  measurement  is  the  least  disturbed  by  topographic  scattering  we  may  expect 
that  at  least  40%  of  the  network  magnitude  variance  shared  by  log(Pmax)  and  log(Pa)  has  little 
or  nothing  to  do  with  topographic  scattering.  That  is  to  say  that  we  can  attribute  at  most  60% 
of  the  log(Pmax)  variance  to  topographic  scattering  and  that  at  least  40%  of  the  variance  must 
be  due  to  other  perturbations.  This  would  place  an  upper  bound  on  the  rms  magnitude  varia¬ 
tion  due  to  topographic  scattering  of  about  0.3  magnitude  units. 
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Additional  evidence  comes  from  the  correlation  of  event  magnitude  residual  patterns  and 
the  "average"  residual  pattern.  The  m^max)  average  residual  pattern  for  the  Ahaggar  test  site 
is  shown  in  Figure  10.  This  residual  pattern  has  an  rms  variation  of  0.30  magnitude  units  and 
accounts  for  40%  of  the  variance  of  mb(Pa),  tr^Ph),  and  m^Pmax)  for  any  of  the  French 
Sahara  shots.  This  would  place  an  upper  limit  of  0.3  rms  magnitude  units  on  the  topographic 
effect  since  the  events  were  located  at  different  places  within  the  mountain.  The  residual  pat¬ 
tern  of  Figure  10  is  consistent  with  focusing-defocusing  of  seismic  energy  across  the  network. 
Clusters  of  large  positive  and  negative  residuals  can  seen  in  the  network  with  rapid  variation  in 
takeoff  angle. 

The  station  residuals  for  both  m^Pa)  and  mb(Pb)  are  50%  correlated  between  RUBIS  and 
SAPHIR  indicating  an  upper  limit  of  0.2  magnitude  units  for  the  topographic  scattering  effect. 
Since  RUBIS  and  SAPHIR  are  expected  to  show  uncorrelated  variations  due  to  topographic 
scattering,  we  can  place  an  upper  bound  of  0.2  magnitude  units  on  the  rms  magnitude  variation 
due  to  topographic  scattering.  Finally,  the  mtCPmax)  station  residuals  are  80%  correlated 
between  SAPHIR  and  RUBIS  and  this  would  suggest  an  upper  bound  of  0.15  rms  magnitude 
units  due  to  topographic  scattering. 

We  are  left  with  a  maximum  log(Pmax)  rms  variation  of  0.15  magnitude  units  across  the 
WWSSN  network  due  to  topographic  scattering,  This  is  in  line  with  what  could  be  expected 
due  to  the  vagaries  of  P+pP  interference  variations  due  to  topographic  scattering.  However,  it 
is  clear  that  the  dominant  contribution  to  the  mb  variance  is  due  to  event-station  perturbations 
that  are  common  between  events  for  the  test  site  and  therefore  are  either  deep  seated  beneath 
the  test  site,  or  are  associated  with  the  stations. 

CONCLUSIONS 
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We  have  made  qualitative  comparisons  of  2-D  LFD  calculations  for  line  sources  under  a 
topographic  feature  with  observations  of  explosions  at  the  French  southern  Sahara  Test  Site. 
The  2-D  simulations  suggest  that  the  linear  elastic  free  surface  (pP)  reflections  could  be 
strongly  affected  by  topography  for  RUBIS  and  SAPH1R.  The  presence  or  absence  of  a  pP 
may  produce  a  0.3  magnitude  variation  in  the  network  mb.  The  topography  could  also  produce 
variations  of  0.3  magnitude  units  in  the  P-wave  teleseismic  magnitudes  by  scattering  of 
Rayleigh-to-P  at  the  free  surface.  However,  since  the  2-D  calculations  probably  over  estimate 
the  P  coda,  this  should  be  considered  an  upper  bound  on  the  magnitude  variations  that  would 
be  introduced  by  such  topography.  Near-source  topographic  scattering  is  expected  to  be  a 
strong  function  of  azimuth  and  takeoff  angle.  These  variations  of  0.3  magnitude  units  would 
be  manifested  in  a  plus  or  minus  0.15  magnitude  variation  about  an  average  and  therefore  be 
responsible  for  at  most  a  0.15  magnitude  variation  about  the  mean. 

2-D  numerical  experiments  suggest  that  spectral  ratios  may  be  diagnostic  of  scattering 
from  topography,  however  simple  spectral  ratio  measures  do  not  show  convincing  evidence  for 
scattering  from  the  topography  at  Taourirt  Tan  Afella  Massif.  This  could  be  due  to  limitations 
of  the  2-D  simulations  or  it  could  be  that  the  bandwidths  of  the  teleseismic  P  waves  are  too 
narrow. 

2-D  numerical  experiments  suggest  that  explosion  line  source  generated  pP  may  be 
focused  or  defocused  by  topography  at  Taourirt  Tan  Afella  Massif.  The  deconvolved  displace¬ 
ment  P  waveforms  at  the  arrays  EKA  and  YKA  and  across  the  network  of  LRSM  stations  do 
show  evidence  for  significant  variation  in  the  free  surface  interaction  with  takeoff  angle  and 
azimuth.  No  simple  pP  is  reliably  observed  and  several  long-petiod  positive  pulses  are 
observed  following  the  the  initial  P  wave  by  1  to  2  seconds.  These  positive  pulses  may  be 
scattered  Rayleigh-to-P  contributions  to  the  P  coda.  The  pP  may  be  defocused  at  some 
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azimuths  and  takeoff  angles  by  the  mountain  topography.  However,  non-linear  interaction  with 
the  free  surface,  such  as  spall,  may  also  be  responsible  for  these  observations. 

Magnitude  statistics  for  the  French  Sahara  explosions  demonstrate  that  there  is  consider¬ 
able  variation  in  the  P  waveforms  observed  at  WWSSN  stations.  The  strong  correlations 
between  log(Pa),  log(Pb)  and  log(Pmax)  and  the  near  source  contribution  to  the  variance  of  the 
logarithmic  ratios  of  log(Pmax/Pa),  log(Pmax/Pb),  and  log(Pb/Pa)  indicate  that  the  contribution 
of  path-receiver  variation  to  the  magnitude  scatter  is  nearly  80%  of  the  magnitude  variance 


regardless  of  the  portion  of  the  P  waveform  that  is  measured.  An  upper  bound  of  0.15  rms 
magnitude  variation  may  be  attributed  to  the  topographic  scattering  from  the  magnitude  statis¬ 
tics.  This  is  in  rough  agreement  with  what  may  be  expected  from  focusing-defocusing  of  pP 
by  topography  or  the  Rayleigh-to-P  scattering  by  the  topography  at  Taourirt  Tan  Afella  Massif. 
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FIGURE  CAPTIONS 

FIGURE  1.  Topographic  map  Taourirt  Tan  Afella  Massif  from  Duclaux  and  Michaud  (1970) 
with  event  locations  from  Faure  (1972)  for  SAPHIR  (S),  RUBIS  (R),  EMERAUDE  (E),  and 
GRENAT  (G).  Contours  are  100  meters.  The  dashed  line  is  the  outcrop  of  the  Taourirt  Tan 
Afella  Massif  granite.  Azimuths  to  several  stations  and  arrays  are  indicated. 

FIGURE  2A.  (ABOVE)  Topographic  profile  with  line  source  locations  indicated  by  triangles. 
FIGURE  2B.  (BELOW)  Logarithmic  amplitude  variations  for  "a",  "b",  and  "max"  P  waves 
radiated  by  the  line  sources  shown  above.  A  0.3  extreme  magnitude  variation  is  produced  in 
the  "max"  phase  due  to  constructive  and/or  destructive  interference  of  the  P+pP  phase. 

FIGURE  3.  Waveforms  for  the  sources  indicated  in  2A.  Left  to  right  and  top  to  bottom 
correspond  each  waveform  to  the  grid  of  source  locations  indicated  in  2A.  A  0.3  maximum 
magnitude  variation  in  the  "max"  P  phase  is  predicted  from  left  to  right  corresponding  to  a 
constant  elevation. 

FIGURE  4A.  North-to-south  topographic  profile  for  Taourirt  Tan  Afella  with  an  incident  20° 
plane  dilatational  wave,  time  -  0.0  sec.  Line  source  locations  are  indicated  by  "S"  and  "R". 
4B,  C,  D,  E  and  F  show  the  amplitude  wave  field  at  time  -  0.125,  0.25,  0.375,  0.5  and  0.675 
sec. 


FIGURE  5A.  Synthetic  far-field  P-wave  seismograms  for  takeoff  angles  of  0,  5,  10,  15,  and  20 
degrees  for  the  line  source  location  indicated  by  "S”  in  Figure  4A.  100  KT  source  with  t* 
=0.45  sec,  and  LRSM  instrument  response. 

FIGURE  5B.  Synthetic  far-field  P-wave  seismograms  for  takeoff  angles  of  0,  5,  10,  15,  and  20 
degrees  for  the  line  source  location  indicated  by  "R"  in  Figure  4A.  100  KT  source  with  t* 
=0.45  sec,  and  LRSM  instrument  response. 

FIGURE  6.  Spectral  ratios  of  RUBIS-to-SAPHIR  at  YKA,  HN-ME,  BL-WV,  LC-NM  and 
EKA.  The  spectral  ratios  are  on  average  flat  with  respect  to  frequency  but  show  variation 
between  stations.  Spectra  were  corrected  for  noise,  and  only  spectral  estimates  with  a  signal- 
to-noise  power  ratio  greater  than  2  are  plotted. 

FIGURE  7.  Equivalent  seismic  sources  from  the  multichannel  deconvolution  of  SAPHIR, 
RUBIS,  GRENAT,  and  EMERAUDE  data  at  the  arrays  EKA  and  YKA  using  the  method  of 
Shumway  and  Der  (1985).  RUBIS  was  poorly  recorded  at  YKA.  Note  the  variation  between 
EKA  and  YKA  for  the  events  EMERAUDE,  GRENAT,  and  SAPHIR.  The  variation  is  the 
greatest  for  SAPHIR  where  a  negative  pulse  that  follows  the  positive  pulse  at  EKA  is  absent  at 
YKA.  The  area  under  the  initial  causal  P-wave  arrival  is  shaded. 


FIGURE  8A.  Deconvolved  seismic  traces  for  RUBIS  at  several  LRSM  stations  arranged  in 
order  of  increasing  azimuth  from  top  to  bottom.  The  area  of  the  initial  causal  P  waveform  is 
indicated  in  each  case.  Secondary  phases  following  the  initial  P  wave  can  be  correlated  over 
several  stations  with  similar  azimuth,  labeled  A  and  (A). 
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FIGURE  8B.  Deconvolved  seismic  traces  for  SAPHIR  at  several  LRSM  stations  arranged  in 
order  of  increasing  azimuth  from  top  to  bottom.  The  area  of  the  initial  P  waveform  is  indi¬ 
cated  in  each  case. 

FIGURE  9A.  Synthetic  teleseismic  P  waves  (takeoff  angle  of  20  degrees)  for  the  2-D 
"SAPHIR"  models  at  azimuths  of  310,  344,  and  0  degrees.  Synthetics  have  only  been  con¬ 
volved  with  an  explosion  source  time  function.  No  distinct  well  defined  elastic  pP  is  apparent 
although,  several  long-period  complicated  negative  pulses  can  be  seen  to  follow  the  initial  P 
wave.  Large  positive  secondary  arrivals  can  be  seen  0.5  sec  following  the  initial  P  wave  at  0 
degrees  azimuth. 

FIGURE  9B.  Synthetic  teleseismic  P  waves  (takeoff  angle  of  20  degrees)  for  the  2-D  "RUBIS" 
models  at  azimuths  of  310,  344,  and  0  degrees.  Synthetics  have  only  been  convolved  with  an 
explosion  source  time  function.  A  distinct  well  defined  elastic  pP  is  only  apparent  for  the 
azimuth  of  0  degrees.  Large  positive  pulses  can  be  seen  about  1  second  after  the  P  wave  for 
the  azimuth  of  344  degrees. 

FIGURE  10.  French  southern  Sahara  test  site  mb  residuals  across  the  WWSSN  network  at  epi- 
central  distances  between  20  and  95  degrees.  Equidistant  polar  projection  with  coast  lines 
(TOP),  and  without  coast  line  (BOTTOM).  Positive  residuals  are  octagons,  negative  residuals 
are  triangles.  The  size  of  the  symbol  is  proportional  to  the  absolute  value  of  the  station  resi¬ 
dual. 
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